Introduction

This report describes the results of a preregistered study available at: https://osf.io/gkd8s.

Note also that this data has been cleaned beforehand. Several datasets over three measurement times were merged (joined) through an inner join so as to keep only participants who at least participated at each step of the study. Missing data will be imputed later on. Duplicates were addressed with the rempsyc::best_duplicate function, which keeps the duplicate with the least amount of missing values, and in case of ties, takes the first occurrence. However, for duplicate participation in the activities and exercises, rather than the first occurrence, the occurrence with the higher completion percentage (% of total activity time) was taken instead.

Packages & Data

…

Packages

library(rempsyc)
library(dplyr)
library(interactions)
library(performance)
library(see)
library(report)
library(datawizard)
library(modelbased)
library(ggplot2)
library(bestNormalize)
library(psych)
library(visdat)
library(missForest)
library(doParallel)
library(ggplot2)
library(emmeans)
library(sjlabelled)
library(eefAnalytics)
library(tidyr)

Data

# Read data
data <- readRDS("Data/finaldataset_n381.rds")

report_participants(data, threshold = 1) %>% cat

381 participants (Mean age = 24.9, SD = 4.1, range: [18, 35], 43.0% missing; Gender: 45.1% women, 11.8% men, 0.00% non-binary, 43.04% missing)

# Allocation ratio
report(data$T1_Group)

x: 3 levels, namely Meditation (n = 65, 17.06%), Reflection (n = 56, 14.70%), Waitlist (n = 96, 25.20%) and missing (n = 164, 43.04%)

report(data$T2_Condition)

x: 2 levels, namely Control (n = 119, 31.23%), Depleted (n = 98, 25.72%) and missing (n = 164, 43.04%)

Data Cleaning

In this section, we are preparing the data for analysis: (a) taking care of preliminary exclusions, (b) checking for and exploring missing values, (d) imputing missing data with missForest, (e) computing scale means, and (f) extracting reliability indices for our scales.

Preliminary exclusions

Here, we report on participant exclusions and corresponding changes on the sample.

…

No Group

First, we want to exclude those who did not get assigned a group. How many people would we exclude?

data %>% 
  filter(is.na(T1_Group)) %>% 
  nrow()
## [1] 164

164 people. Investigation of these data suggests these participants did not complete the questionnaires and did not reach the group assignment. It is safe to exclude them, so let’s do so.

data <- data %>% 
  filter(!is.na(T1_Group))

report_participants(data, threshold = 1) %>% cat
## 217 participants (Mean age = 24.9, SD = 4.1, range: [18, 35]; Gender: 79.3% women, 20.7% men, 0.00% non-binary)
# Allocation ratio
report(data$T1_Group)
## x: 3 levels, namely Meditation (n = 65, 29.95%), Reflection (n = 56, 25.81%)
## and Waitlist (n = 96, 44.24%)
report(data$T2_Condition)
## x: 2 levels, namely Control (n = 119, 54.84%) and Depleted (n = 98, 45.16%)

Participation %

Second, we know that we only want to keep participants who had a participation level of at least 2/3 of all activities and exercises. Let’s see the distribution of participants’ participation, by group. However, we do not want to exclude participants from the control group, so we will give them an artificial 100% participation rate. How many people would we exclude?

data <- data %>% 
  mutate(part.percent = convert_na_to(part.percent, 1))

data %>% 
  filter(part.percent < 2/3) %>% 
  count(T1_Group)
T1_Group n
Meditation 6
Reflection 1

7 people. Let’s exclude them.

data2 <- data

data <- data %>% 
  filter(part.percent >= 2/3)

report_participants(data, threshold = 1) %>% cat
## 210 participants (Mean age = 24.8, SD = 4.0, range: [18, 35]; Gender: 80.5% women, 19.5% men, 0.00% non-binary)
# Allocation ratio
report(data$T1_Group)
## x: 3 levels, namely Meditation (n = 59, 28.10%), Reflection (n = 55, 26.19%)
## and Waitlist (n = 96, 45.71%)
report(data$T2_Condition)
## x: 2 levels, namely Control (n = 113, 53.81%) and Depleted (n = 97, 46.19%)

Attention Checks

Let’s also exclude those who failed 2 or more attention checks (i.e., keep with those with a score of two or more).

data <- data %>% 
    mutate(att_check = rowSums(
      select(., T1_attention1, T2_attention1, T3_attention1), na.rm = TRUE))

data %>% 
  count(att_check)
att_check n
0 3
1 4
2 13
3 190

There’s 7 more exclusions here.

data <- data %>% 
  filter(att_check >= 2)

Demographics

Here, we report various sample demographics.

…

Sample

report_participants(data, threshold = 1) %>% cat

203 participants (Mean age = 24.8, SD = 3.9, range: [18, 35]; Gender: 81.3% women, 18.7% men, 0.00% non-binary)

report_participants(data, threshold = 1, group = "T1_Group") %>% cat

For the ‘T1_Group - Meditation’ group: 58 participants (Mean age = 25.2, SD = 4.5, range: [18, 35]), for the ‘T1_Group - Reflection’ group: 53 participants (Mean age = 25.1, SD = 3.7, range: [19, 35]) and for the ‘T1_Group - Waitlist’ group: 92 participants (Mean age = 24.4, SD = 3.7, range: [18, 35])

# Allocation ratio
report(data$T1_Group)

x: 3 levels, namely Meditation (n = 58, 28.57%), Reflection (n = 53, 26.11%) and Waitlist (n = 92, 45.32%)

report(data$T2_Condition)

x: 2 levels, namely Control (n = 113, 55.67%) and Depleted (n = 90, 44.33%)

Recruitment

“How did you hear about the study?”

get_label(data$T1_recruitment) %>% cat
## Comment avez-vous entendu parler de l'étude? - Selected Choice
report(data$T1_recruitment)
## x: 8 levels, namely Email (n = 49, 24.14%), Facebook (n = 116, 57.14%), google
## (n = 5, 2.46%), Kijiji (n = 8, 3.94%), Other (n = 0, 0.00%), SQRP bulletin
## board (n = 10, 4.93%), SSMU Marketplace (n = 6, 2.96%) and Word of mouth (n =
## 9, 4.43%)
data %>% 
    count(T1_recruitment, sort = TRUE)
T1_recruitment n
Facebook 116
Email 49
SQRP bulletin board 10
Word of mouth 9
Kijiji 8
SSMU Marketplace 6
google 5

Age

data %>% 
  nice_density("age", histogram = TRUE)

Gender

data %>% 
    count(gender, sort = TRUE)
gender n
Female 165
Male 38

Psychology Classes

“Have you already completed a psychology course?”

get_label(data$T1_psycho.class) %>% cat
## Avez-vous déjà suivi un cours de psychologie?
data %>% 
    count(T1_psycho.class, sort = TRUE)
T1_psycho.class n
Yes 124
No 79

Virtual Reality

“Have you ever tried virtual reality?”

get_label(data$T1_virtual.reality) %>% cat
## Avez-vous déjà essayé une sorte de réalité virtuelle?
data %>% 
    count(T1_virtual.reality, sort = TRUE)
T1_virtual.reality n
No 137
Yes 66

Meditation Experience

“What is your meditation experience?”

get_label(data$T1_medi.experience) %>% cat
## Quelle expérience de méditation avez-vous?
data %>% 
    count(T1_medi.experience, sort = TRUE, .drop = FALSE)
T1_medi.experience n
< 5 hours 161
Between 5 and 10 hours 42
> 10 hours 0

Disorders

“Do you have a history of psychiatric or neurological disorders?”

get_label(data$T1_disorders) %>% cat
## Avez-vous des antécédents de troubles psychiatriques ou neurologiques?
data %>% 
    count(T1_disorders, sort = TRUE)
T1_disorders n
No 203

Vision

“Do you have normal or corrected-to-normal vision?”

get_label(data$T1_vision) %>% cat
## Avez-vous une vision normale ou corrigée à la normale?
data %>% 
    count(T1_vision, sort = TRUE)
T1_vision n
Yes 203

Phone

“Do you own a smart phone of type iPhone or Android?”

get_label(data$T1_phone) %>% cat
## Possédez-vous un téléphone intelligent de type iPhone ou Android?
data %>% 
    count(T1_phone, sort = TRUE)
T1_phone n
Yes 201
No 2

Live in Quebec

“Do you live in Quebec at this time?”

get_label(data$T1_quebec) %>% cat
## Habitez-vous au Québec présentement?
data %>% 
    count(T1_quebec, sort = TRUE)
T1_quebec n
NA 168
Yes 34
No 1

Student

get_label(data$T1_student) %>% cat
## Êtes-vous étudiant.e?
data %>% 
    count(T1_student, sort = TRUE)
T1_student n
Student 160
Non-student 43 
get_label(data$T1_student.program_cat) %>% cat
## Dans quel programme étudiez-vous?
report(data$T1_student.program)
## x: 111 entries, such as psychologie (11.33%); baccalauréat en psychologie
## (2.96%); criminologie (2.96%) and 108 others (43 missing)
data %>% 
    count(T1_student.program_cat, sort = TRUE) %>% 
  filter(n > 3)
T1_student.program_cat n
Psychology 45
NA 43
Other 29
Biology 7
Teaching 7
Biomedical Sciences 6
Law 5
Nursing 5
Social Work 5
Arts 4
History 4
Neuroscience 4
Sexology 4

Workplace

get_label(data$T1_workplace) %>% cat
## Dans quel domaine travaillez-vous?
report(data$T1_workplace)
## x: 39 entries, such as informatique (1.48%); administration (0.99%);
## psychologie (0.99%) and 36 others (160 missing)
data %>% 
    count(T1_workplace_cat, sort = TRUE) %>% 
  filter(n > 1)
T1_workplace_cat n
NA 160
Other 12
Health 5
Information and Technology 5
Psychology 5
Admin 4
Arts & Culture 4

Meditation Practice

Participants were asked questions about their meditation practice during the 13 weeks.

get_label(data$T3_post.medipractice) %>% cat
## Avez-vous pratiqué la méditation au cours des 13 dernières semaines?
data %>% 
    count(T3_post.medipractice, sort = TRUE)
T3_post.medipractice n
Non 107
Oui 95
NA 1 
get_label(data$T3_medipractice.which) %>% cat
## Quel type de méditation avez-vous pratiqué? - Selected Choice
report(data$T3_medipractice.which)
## x: 21 entries, such as NA (53.20%); mindfulness (20.69%); LKM (7.88%) and 18
## others (0 missing)
data %>% 
    count(T3_medipractice.which, sort = TRUE)
T3_medipractice.which n
NA 108
mindfulness 42
LKM 16
mindfulness, LKM 16
chakras 3
mindfulness, LKM, vipassana 2
other, (Guidée) 2
LKM, chakras 1
LKM, chakras, transcendental 1
mindfulness, LKM, chakras 1
mindfulness, LKM, chakras, other, (amour-propre et healing) 1
mindfulness, LKM, other, (J’ai fait les 42 jours de méditation par l’étude et de la méditation seul (avec bol)) 1
mindfulness, LKM, other, (en marchant) 1
mindfulness, chakras 1
mindfulness, other, (Respiration) 1
mindfulness, transcendental 1
other, (Générique, juste se coucher et se vider l’esprit) 1
other, (Libre pensée / se-laisser-porter) 1
other, (Relaxation) 1
other, (Respiration-visualisation) 1
vipassana 1 
get_label(data$T3_medipractice.time) %>% cat
## À quelle fréquence avez-vous médité par semaine?
data %>% 
    count(T3_medipractice.time, sort = TRUE)
T3_medipractice.time n
NA 108
10-20 min 35
< 10 min 33
20-30 min 10
1-2 h 8
30-60 min 8
> 2 h 1 
get_label(data$T3_choice.medicomp) %>% cat
## Aimeriez-vous recevoir le programme de méditation gratuitement?
data %>% 
    count(T3_choice.medicomp, sort = TRUE)
T3_choice.medicomp n
Oui 104
Non 93
already received 5
NA 1

Explore missing data

Here, we explore missing data, first by item and questionnaire, then using the visdat package, and finally using Little’s MCAR test.

…

Missing items

# Check for nice_na
nice_na(data, scales = c(
  "T1_BSCS", "T1_BAQ", "T1_NOBAGS", "T1_attitude", "T1_dehumanization", 
  "T1_WHS", "T1_CLS", "T2_NOBAGS", "T2_attitude", "T2_dehumanization", 
  "T2_SMS5", "T2_PANAS", "T2_WHS", "T2_CLS", "T2_charity", "T3_NOBAGS", 
  "T3_attitude", "T3_dehumanization", "T3_WHS", "T3_CLS"))
var items na cells na_percent na_max na_max_percent all_na
T1_BSCS_1:T1_BSCS_7 7 0 1421 0.00 0 0.00 0
T1_BAQ_1:T1_BAQ_13 13 1 2639 0.04 1 7.69 0
T1_NOBAGS.1_1:T1_NOBAGS.16_1 20 0 4060 0.00 0 0.00 0
T1_attitude_1:T1_attitude_9 9 0 1827 0.00 0 0.00 0
T1_dehumanization_1:T1_dehumanization_9 9 1 1827 0.05 1 11.11 0
T1_WHS_1:T1_WHS_6 6 0 1218 0.00 0 0.00 0
T1_CLS_1:T1_CLS_21 21 0 4263 0.00 0 0.00 0
T2_NOBAGS.1_1:T2_NOBAGS.16_1 20 0 4060 0.00 0 0.00 0
T2_attitude_1:T2_attitude_9 9 0 1827 0.00 0 0.00 0
T2_dehumanization_1:T2_dehumanization_9 9 1 1827 0.05 1 11.11 0
T2_SMS5_1:T2_SMS5_6 6 0 1218 0.00 0 0.00 0
T2_PANAS_1:T2_PANAS_10 10 0 2030 0.00 0 0.00 0
T2_WHS_1:T2_WHS_6 6 0 1218 0.00 0 0.00 0
T2_CLS_1:T2_CLS_21 21 0 4263 0.00 0 0.00 0
T2_charity.moisson1_1:T2_charity.armee2_1 48 62 9744 0.64 22 45.83 0
T3_NOBAGS.1_1:T3_NOBAGS.16_1 20 0 4060 0.00 0 0.00 0
T3_attitude_1:T3_attitude_9 9 0 1827 0.00 0 0.00 0
T3_dehumanization_1:T3_dehumanization_9 9 0 1827 0.00 0 0.00 0
T3_WHS_1:T3_WHS_6 6 0 1218 0.00 0 0.00 0
T3_CLS_1:T3_CLS_21 21 21 4263 0.49 21 100.00 1
Total 389 6763 78967 8.56 89 22.88 0

A few items are missing here and there.

Patterns of missing data

Let’s check for patterns of missing data.

# Smaller subset of data for easier inspection
data %>%
  # select(manualworkerId:att_check2_raw, 
  #        condition:condition_dum) %>%
  vis_miss

Little’s MCAR test

# Let's use Little's MCAR test to confirm
# We have to proceed by "scale" because the function can only
# support 30 variables max at a time
library(naniar)
## Warning: package 'naniar' was built under R version 4.2.3
# We only check for the variable that had missing data, charity

# Have to divide this one in two because it is too large for the function
data %>% 
  select(T2_charity.moisson1_1:T2_charity.suzuki2_1) %>% 
  mcar_test
statistic df p.value missing.patterns
214.1051 181 0.0466264 8 
data %>% 
  select(T2_charity.conserv1_1:T2_charity.armee2_1) %>% 
  mcar_test
statistic df p.value missing.patterns
69.62351 61 0.2099961 6

Impute missing data

Here, we impute missing data with the missForest package, as it is one of the best imputation methods.

…

Imputation

Details

Why impute the data? van Ginkel explains,

Regardless of the missingness mechanism, multiple imputation is always to be preferred over listwise deletion. Under MCAR it is preferred because it results in more statistical power, under MAR it is preferred because besides more power it will give unbiased results whereas listwise deletion may not, and under NMAR it is also the preferred method because it will give less biased results than listwise deletion.

van Ginkel, J. R., Linting, M., Rippe, R. C. A., & van der Voort, A. (2020). Rebutting existing misconceptions about multiple imputation as a method for handling missing data. Journal of Personality Assessment, 102(3), 297-308. https://doi.org/10.1080/00223891.2018.1530680

Why missForest? It outperforms other imputation methods, including the popular MICE (multiple imputation by chained equations). You also don’t end up with several datasets, which makes it easier for following analyses. Finally, it can be applied to mixed data types (missings in numeric & categorical variables).

Waljee, A. K., Mukherjee, A., Singal, A. G., Zhang, Y., Warren, J., Balis, U., … & Higgins, P. D. (2013). Comparison of imputation methods for missing laboratory data in medicine. BMJ open, 3(8), e002847. https://doi.org/10.1093/bioinformatics/btr597

Stekhoven, D. J., & Bühlmann, P. (2012). MissForest—non-parametric missing value imputation for mixed-type data. Bioinformatics, 28(1), 112-118. https://doi.org/10.1093/bioinformatics/btr597

Scale Means & Reliability

Now that we have imputed the missing data, we are ready to calculate our scale means. After reversing our items, we can then get the alphas and omegas for our different scales.

…

Self-control

# Reverse code items 2, 4, 6, 7
data <- data %>% 
  mutate(across(contains("BSCS"), .names = "{col}r"))

data <- data %>% 
  mutate(across(ends_with(paste("BSCS", c(2, 4, 6, 7), sep = "_")), 
                ~nice_reverse(.x, 5), .names = "{col}r"))

# Get mean BSCS
data <- data %>% 
  mutate(T1_BSCS = rowMeans(pick(T1_BSCS_1r:T1_BSCS_7r)))

# Get alpha & omega
data %>% 
  select(T1_BSCS_1r:T1_BSCS_7r) %>% 
  omega(nfactors = 1)
## Loading required namespace: GPArotation
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.75 
## G.6:                   0.75 
## Omega Hierarchical:    0.75 
## Omega H asymptotic:    1 
## Omega Total            0.75 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##               g  F1*   h2   u2 p2
## T1_BSCS_1r 0.54      0.29 0.71  1
## T1_BSCS_2r 0.55      0.30 0.70  1
## T1_BSCS_3r 0.69      0.47 0.53  1
## T1_BSCS_4r 0.59      0.35 0.65  1
## T1_BSCS_5r 0.48      0.23 0.77  1
## T1_BSCS_6r 0.61      0.37 0.63  1
## T1_BSCS_7r 0.37      0.14 0.86  1
## 
## With Sums of squares  of:
##   g F1* 
## 2.1 0.0 
## 
## general/max  25873775754603412   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 14  and the fit is  0.39 
## The number of observations was  203  with Chi Square =  77.58  with prob <  0.000000000079
## The root mean square of the residuals is  0.09 
## The df corrected root mean square of the residuals is  0.12
## RMSEA index =  0.149  and the 10 % confidence intervals are  0.118 0.183
## BIC =  3.19
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 14  and the fit is  0.39 
## The number of observations was  203  with Chi Square =  77.58  with prob <  0.000000000079
## The root mean square of the residuals is  0.09 
## The df corrected root mean square of the residuals is  0.12 
## 
## RMSEA index =  0.149  and the 10 % confidence intervals are  0.118 0.183
## BIC =  3.19 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.88   0
## Multiple R square of scores with factors      0.77   0
## Minimum correlation of factor score estimates 0.54  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.75 0.75
## Omega general for total scores and subscales  0.75 0.75
## Omega group for total scores and subscales    0.00 0.00

Trait aggression

# Reverse code item 7
data <- data %>% 
  mutate(across(contains("BAQ"), .names = "{col}r"))

data <- data %>% 
  mutate(across(T1_BAQ_7, ~nice_reverse(.x, 7), .names = "{col}r"))

# Get mean BAQ
data <- data %>% 
  mutate(T1_BAQ = rowMeans(pick(T1_BAQ_1r:T1_BAQ_12r)))

# Get alpha & omega
data %>% 
  select(T1_BAQ_1r:T1_BAQ_12r) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.78 
## G.6:                   0.82 
## Omega Hierarchical:    0.79 
## Omega H asymptotic:    1 
## Omega Total            0.79 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##               g  F1*   h2   u2 p2
## T1_BAQ_1r  0.65      0.42 0.58  1
## T1_BAQ_2r  0.57      0.32 0.68  1
## T1_BAQ_3r  0.56      0.31 0.69  1
## T1_BAQ_4r            0.03 0.97  1
## T1_BAQ_5r  0.56      0.31 0.69  1
## T1_BAQ_6r  0.46      0.21 0.79  1
## T1_BAQ_7r  0.41      0.17 0.83  1
## T1_BAQ_8r  0.62      0.39 0.61  1
## T1_BAQ_9r  0.72      0.52 0.48  1
## T1_BAQ_10r 0.33      0.11 0.89  1
## T1_BAQ_11r 0.45      0.20 0.80  1
## T1_BAQ_12r 0.28      0.08 0.92  1
## 
## With Sums of squares  of:
##   g F1* 
## 3.1 0.0 
## 
## general/max  36926570500479704   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 54  and the fit is  1.38 
## The number of observations was  203  with Chi Square =  270.72  with prob <  0.0000000000000000000000000000013
## The root mean square of the residuals is  0.12 
## The df corrected root mean square of the residuals is  0.13
## RMSEA index =  0.141  and the 10 % confidence intervals are  0.125 0.158
## BIC =  -16.19
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 54  and the fit is  1.38 
## The number of observations was  203  with Chi Square =  270.72  with prob <  0.0000000000000000000000000000013
## The root mean square of the residuals is  0.12 
## The df corrected root mean square of the residuals is  0.13 
## 
## RMSEA index =  0.141  and the 10 % confidence intervals are  0.125 0.158
## BIC =  -16.19 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.91   0
## Multiple R square of scores with factors      0.83   0
## Minimum correlation of factor score estimates 0.66  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.79 0.79
## Omega general for total scores and subscales  0.79 0.79
## Omega group for total scores and subscales    0.00 0.00

Explicit Attitudes

data <- data %>% 
  mutate(T1_attitude = rowMeans(pick(T1_attitude_1:T1_attitude_9)),
         T2_attitude = rowMeans(pick(T2_attitude_1:T2_attitude_9)),
         T3_attitude = rowMeans(pick(T3_attitude_1:T3_attitude_9)))

# Get alpha & omega
data %>% 
  select(T1_attitude_1:T1_attitude_9) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in cov2cor(t(w) %*% r %*% w): diag(.) had 0 or NA entries; non-finite
## result is doubtful
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.93 
## G.6:                   0.94 
## Omega Hierarchical:    0.94 
## Omega H asymptotic:    1 
## Omega Total            0.94 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                  g  F1*   h2   u2 p2
## T1_attitude_1 0.92      0.84 0.16  1
## T1_attitude_2 0.78      0.61 0.39  1
## T1_attitude_3 0.91      0.82 0.18  1
## T1_attitude_4 0.85      0.73 0.27  1
## T1_attitude_5 0.83      0.68 0.32  1
## T1_attitude_6 0.83      0.68 0.32  1
## T1_attitude_7 0.39      0.15 0.85  1
## T1_attitude_8 0.83      0.69 0.31  1
## T1_attitude_9 0.69      0.47 0.53  1
## 
## With Sums of squares  of:
##   g F1* 
## 5.7 0.0 
## 
## general/max  Inf   max/min =   NaN
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 27  and the fit is  0.73 
## The number of observations was  203  with Chi Square =  143.81  with prob <  0.0000000000000000068
## The root mean square of the residuals is  0.05 
## The df corrected root mean square of the residuals is  0.06
## RMSEA index =  0.146  and the 10 % confidence intervals are  0.123 0.17
## BIC =  0.35
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27  and the fit is  0.73 
## The number of observations was  203  with Chi Square =  143.81  with prob <  0.0000000000000000068
## The root mean square of the residuals is  0.05 
## The df corrected root mean square of the residuals is  0.06 
## 
## RMSEA index =  0.146  and the 10 % confidence intervals are  0.123 0.17
## BIC =  0.35 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.98   0
## Multiple R square of scores with factors      0.96   0
## Minimum correlation of factor score estimates 0.91  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.94 0.94
## Omega general for total scores and subscales  0.94 0.94
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T2_attitude_1:T2_attitude_9) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor

## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## diag(.) had 0 or NA entries; non-finite result is doubtful
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.96 
## G.6:                   0.97 
## Omega Hierarchical:    0.96 
## Omega H asymptotic:    1 
## Omega Total            0.96 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                  g  F1*   h2   u2 p2
## T2_attitude_1 0.92      0.85 0.15  1
## T2_attitude_2 0.87      0.75 0.25  1
## T2_attitude_3 0.94      0.89 0.11  1
## T2_attitude_4 0.88      0.77 0.23  1
## T2_attitude_5 0.91      0.82 0.18  1
## T2_attitude_6 0.91      0.83 0.17  1
## T2_attitude_7 0.56      0.32 0.68  1
## T2_attitude_8 0.93      0.86 0.14  1
## T2_attitude_9 0.84      0.71 0.29  1
## 
## With Sums of squares  of:
##   g F1* 
## 6.8 0.0 
## 
## general/max  Inf   max/min =   NaN
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 27  and the fit is  0.83 
## The number of observations was  203  with Chi Square =  164.1  with prob <  0.0000000000000000000013
## The root mean square of the residuals is  0.03 
## The df corrected root mean square of the residuals is  0.03
## RMSEA index =  0.158  and the 10 % confidence intervals are  0.136 0.182
## BIC =  20.65
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27  and the fit is  0.83 
## The number of observations was  203  with Chi Square =  164.1  with prob <  0.0000000000000000000013
## The root mean square of the residuals is  0.03 
## The df corrected root mean square of the residuals is  0.03 
## 
## RMSEA index =  0.158  and the 10 % confidence intervals are  0.136 0.182
## BIC =  20.65 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.99   0
## Multiple R square of scores with factors      0.98   0
## Minimum correlation of factor score estimates 0.95  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.96 0.96
## Omega general for total scores and subscales  0.96 0.96
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T3_attitude_1:T3_attitude_9) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.96 
## G.6:                   0.97 
## Omega Hierarchical:    0.96 
## Omega H asymptotic:    1 
## Omega Total            0.96 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                  g  F1*   h2   u2 p2
## T3_attitude_1 0.96      0.93 0.07  1
## T3_attitude_2 0.83      0.68 0.32  1
## T3_attitude_3 0.94      0.88 0.12  1
## T3_attitude_4 0.86      0.74 0.26  1
## T3_attitude_5 0.92      0.85 0.15  1
## T3_attitude_6 0.88      0.78 0.22  1
## T3_attitude_7 0.61      0.37 0.63  1
## T3_attitude_8 0.90      0.81 0.19  1
## T3_attitude_9 0.79      0.62 0.38  1
## 
## With Sums of squares  of:
##   g F1* 
## 6.7 0.0 
## 
## general/max  119952105224445488   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 27  and the fit is  0.98 
## The number of observations was  203  with Chi Square =  193.76  with prob <  0.0000000000000000000000000038
## The root mean square of the residuals is  0.04 
## The df corrected root mean square of the residuals is  0.04
## RMSEA index =  0.174  and the 10 % confidence intervals are  0.152 0.198
## BIC =  50.31
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27  and the fit is  0.98 
## The number of observations was  203  with Chi Square =  193.76  with prob <  0.0000000000000000000000000038
## The root mean square of the residuals is  0.04 
## The df corrected root mean square of the residuals is  0.04 
## 
## RMSEA index =  0.174  and the 10 % confidence intervals are  0.152 0.198
## BIC =  50.31 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.99   0
## Multiple R square of scores with factors      0.98   0
## Minimum correlation of factor score estimates 0.96  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.96 0.96
## Omega general for total scores and subscales  0.96 0.96
## Omega group for total scores and subscales    0.00 0.00

Dehumanization

data <- data %>% 
  mutate(T1_dehumanization = rowMeans(pick(T1_dehumanization_1:T1_dehumanization_9)),
         T2_dehumanization = rowMeans(pick(T2_dehumanization_1:T2_dehumanization_9)),
         T3_dehumanization = rowMeans(pick(T3_dehumanization_1:T3_dehumanization_9)))

# Get alpha & omega
data %>% 
  select(T1_dehumanization_1:T1_dehumanization_9) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.94 
## G.6:                   0.96 
## Omega Hierarchical:    0.96 
## Omega H asymptotic:    1 
## Omega Total            0.96 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                        g  F1*   h2   u2 p2
## T1_dehumanization_1 0.91      0.83 0.17  1
## T1_dehumanization_2 0.93      0.87 0.13  1
## T1_dehumanization_3 0.93      0.86 0.14  1
## T1_dehumanization_4 0.93      0.86 0.14  1
## T1_dehumanization_5 0.95      0.90 0.10  1
## T1_dehumanization_6 0.87      0.76 0.24  1
## T1_dehumanization_7           0.00 1.00  1
## T1_dehumanization_8 0.92      0.85 0.15  1
## T1_dehumanization_9 0.81      0.66 0.34  1
## 
## With Sums of squares  of:
##   g F1* 
## 6.6 0.0 
## 
## general/max  337359029506734720   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 27  and the fit is  1.24 
## The number of observations was  203  with Chi Square =  244.36  with prob <  0.00000000000000000000000000000000000069
## The root mean square of the residuals is  0.03 
## The df corrected root mean square of the residuals is  0.04
## RMSEA index =  0.199  and the 10 % confidence intervals are  0.177 0.223
## BIC =  100.9
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27  and the fit is  1.24 
## The number of observations was  203  with Chi Square =  244.36  with prob <  0.00000000000000000000000000000000000069
## The root mean square of the residuals is  0.03 
## The df corrected root mean square of the residuals is  0.04 
## 
## RMSEA index =  0.199  and the 10 % confidence intervals are  0.177 0.223
## BIC =  100.9 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.99   0
## Multiple R square of scores with factors      0.98   0
## Minimum correlation of factor score estimates 0.96  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.96 0.96
## Omega general for total scores and subscales  0.96 0.96
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T2_dehumanization_1:T2_dehumanization_9) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.96 
## G.6:                   0.97 
## Omega Hierarchical:    0.97 
## Omega H asymptotic:    1 
## Omega Total            0.97 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                        g  F1*   h2   u2 p2
## T2_dehumanization_1 0.98      0.96 0.04  1
## T2_dehumanization_2 0.94      0.89 0.11  1
## T2_dehumanization_3 0.94      0.89 0.11  1
## T2_dehumanization_4 0.92      0.84 0.16  1
## T2_dehumanization_5 0.96      0.92 0.08  1
## T2_dehumanization_6 0.92      0.85 0.15  1
## T2_dehumanization_7           0.00 1.00  1
## T2_dehumanization_8 0.95      0.91 0.09  1
## T2_dehumanization_9 0.94      0.88 0.12  1
## 
## With Sums of squares  of:
##   g F1* 
## 7.1 0.0 
## 
## general/max  391733323604288832   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 27  and the fit is  1.09 
## The number of observations was  203  with Chi Square =  215.86  with prob <  0.00000000000000000000000000000023
## The root mean square of the residuals is  0.02 
## The df corrected root mean square of the residuals is  0.02
## RMSEA index =  0.186  and the 10 % confidence intervals are  0.163 0.21
## BIC =  72.41
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27  and the fit is  1.09 
## The number of observations was  203  with Chi Square =  215.86  with prob <  0.00000000000000000000000000000023
## The root mean square of the residuals is  0.02 
## The df corrected root mean square of the residuals is  0.02 
## 
## RMSEA index =  0.186  and the 10 % confidence intervals are  0.163 0.21
## BIC =  72.41 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.99   0
## Multiple R square of scores with factors      0.99   0
## Minimum correlation of factor score estimates 0.98  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.97 0.97
## Omega general for total scores and subscales  0.97 0.97
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T3_dehumanization_1:T3_dehumanization_9) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in cov2cor(t(w) %*% r %*% w): diag(.) had 0 or NA entries; non-finite
## result is doubtful
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.95 
## G.6:                   0.97 
## Omega Hierarchical:    0.96 
## Omega H asymptotic:    1 
## Omega Total            0.96 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                        g  F1*   h2   u2 p2
## T3_dehumanization_1 0.95      0.90 0.10  1
## T3_dehumanization_2 0.93      0.87 0.13  1
## T3_dehumanization_3 0.93      0.87 0.13  1
## T3_dehumanization_4 0.95      0.91 0.09  1
## T3_dehumanization_5 0.95      0.91 0.09  1
## T3_dehumanization_6 0.89      0.78 0.22  1
## T3_dehumanization_7           0.01 0.99  1
## T3_dehumanization_8 0.93      0.86 0.14  1
## T3_dehumanization_9 0.84      0.70 0.30  1
## 
## With Sums of squares  of:
##   g F1* 
## 6.8 0.0 
## 
## general/max  Inf   max/min =   NaN
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 27  and the fit is  2.32 
## The number of observations was  203  with Chi Square =  457.5  with prob <  0.000000000000000000000000000000000000000000000000000000000000000000000000000000087
## The root mean square of the residuals is  0.03 
## The df corrected root mean square of the residuals is  0.04
## RMSEA index =  0.28  and the 10 % confidence intervals are  0.259 0.304
## BIC =  314.04
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27  and the fit is  2.32 
## The number of observations was  203  with Chi Square =  457.5  with prob <  0.000000000000000000000000000000000000000000000000000000000000000000000000000000087
## The root mean square of the residuals is  0.03 
## The df corrected root mean square of the residuals is  0.04 
## 
## RMSEA index =  0.28  and the 10 % confidence intervals are  0.259 0.304
## BIC =  314.04 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.99   0
## Multiple R square of scores with factors      0.98   0
## Minimum correlation of factor score estimates 0.96  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.96 0.96
## Omega general for total scores and subscales  0.96 0.96
## Omega group for total scores and subscales    0.00 0.00

Aggression Attitude (NOBAGS)

data <- data %>% 
  mutate(across(contains("NOBAGS"), .names = "{col}r"))

# Reverse code NOBAGS (items 1:2, 5:6, 10,12, 14:16, 20)
data <- data %>%
  mutate(across(ends_with(paste0("NOBAGS.", c(
    "1_1", "1_2", "3_1", "3_2", "6_1", "8_1", "10_1", "12_1", "16_1"))), 
    ~nice_reverse(.x, 4), .names = "{col}r"))

# Get mean NOBAGS
data <- data %>% 
  mutate(T1_NOBAGS = rowMeans(pick(T1_NOBAGS.1_1r:T1_NOBAGS.16_1r)),
         T2_NOBAGS = rowMeans(pick(T2_NOBAGS.1_1r:T2_NOBAGS.16_1r)),
         T3_NOBAGS = rowMeans(pick(T3_NOBAGS.1_1r:T3_NOBAGS.16_1r)))

# Get alpha & omega
data %>% 
  select(T1_NOBAGS.1_1r:T1_NOBAGS.16_1r) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.85 
## G.6:                   0.91 
## Omega Hierarchical:    0.84 
## Omega H asymptotic:    0.99 
## Omega Total            0.85 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                     g  F1*   h2   u2 p2
## T1_NOBAGS.1_1r   0.63      0.40 0.60  1
## T1_NOBAGS.1_2r   0.52      0.27 0.73  1
## T1_NOBAGS.2_1r   0.53      0.28 0.72  1
## T1_NOBAGS.2_2r   0.46      0.21 0.79  1
## T1_NOBAGS.3_1r   0.64      0.41 0.59  1
## T1_NOBAGS.3_2r   0.50      0.25 0.75  1
## T1_NOBAGS.4_1r   0.67      0.45 0.55  1
## T1_NOBAGS.4_2r   0.47      0.22 0.78  1
## T1_NOBAGS.5_1r   0.51      0.26 0.74  1
## T1_NOBAGS.6_1r   0.42      0.18 0.82  1
## T1_NOBAGS.7_1r   0.57      0.33 0.67  1
## T1_NOBAGS.8_1r   0.45      0.20 0.80  1
## T1_NOBAGS.9_1r   0.44      0.20 0.80  1
## T1_NOBAGS.10_1r  0.34      0.12 0.88  1
## T1_NOBAGS.11_1r- 0.27      0.07 0.93  1
## T1_NOBAGS.12_1r  0.30      0.09 0.91  1
## T1_NOBAGS.13_1r  0.43      0.18 0.82  1
## T1_NOBAGS.14_1r  0.45      0.21 0.79  1
## T1_NOBAGS.15_1r  0.47      0.22 0.78  1
## T1_NOBAGS.16_1r            0.03 0.97  1
## 
## With Sums of squares  of:
##   g F1* 
## 4.6 0.0 
## 
## general/max  14975690996060718   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 170  and the fit is  6.58 
## The number of observations was  203  with Chi Square =  1276.17  with prob <  0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011
## The root mean square of the residuals is  0.16 
## The df corrected root mean square of the residuals is  0.17
## RMSEA index =  0.179  and the 10 % confidence intervals are  0.17 0.189
## BIC =  372.93
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 170  and the fit is  6.58 
## The number of observations was  203  with Chi Square =  1276.17  with prob <  0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011
## The root mean square of the residuals is  0.16 
## The df corrected root mean square of the residuals is  0.17 
## 
## RMSEA index =  0.179  and the 10 % confidence intervals are  0.17 0.189
## BIC =  372.93 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.93   0
## Multiple R square of scores with factors      0.87   0
## Minimum correlation of factor score estimates 0.73  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.85 0.84
## Omega general for total scores and subscales  0.84 0.84
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T2_NOBAGS.1_1r:T2_NOBAGS.16_1r) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.88 
## G.6:                   0.94 
## Omega Hierarchical:    0.88 
## Omega H asymptotic:    0.99 
## Omega Total            0.88 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                     g  F1*   h2   u2 p2
## T2_NOBAGS.1_1r   0.53      0.28 0.72  1
## T2_NOBAGS.1_2r   0.41      0.17 0.83  1
## T2_NOBAGS.2_1r   0.68      0.46 0.54  1
## T2_NOBAGS.2_2r   0.69      0.47 0.53  1
## T2_NOBAGS.3_1r   0.63      0.40 0.60  1
## T2_NOBAGS.3_2r   0.54      0.29 0.71  1
## T2_NOBAGS.4_1r   0.64      0.41 0.59  1
## T2_NOBAGS.4_2r   0.61      0.37 0.63  1
## T2_NOBAGS.5_1r   0.68      0.47 0.53  1
## T2_NOBAGS.6_1r   0.45      0.20 0.80  1
## T2_NOBAGS.7_1r   0.68      0.46 0.54  1
## T2_NOBAGS.8_1r   0.42      0.18 0.82  1
## T2_NOBAGS.9_1r   0.46      0.21 0.79  1
## T2_NOBAGS.10_1r  0.37      0.13 0.87  1
## T2_NOBAGS.11_1r- 0.38      0.15 0.85  1
## T2_NOBAGS.12_1r  0.46      0.21 0.79  1
## T2_NOBAGS.13_1r  0.47      0.22 0.78  1
## T2_NOBAGS.14_1r  0.47      0.22 0.78  1
## T2_NOBAGS.15_1r  0.42      0.18 0.82  1
## T2_NOBAGS.16_1r  0.42      0.18 0.82  1
## 
## With Sums of squares  of:
##   g F1* 
## 5.7 0.0 
## 
## general/max  17017814006844390   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 170  and the fit is  7.91 
## The number of observations was  203  with Chi Square =  1532.35  with prob <  0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012
## The root mean square of the residuals is  0.18 
## The df corrected root mean square of the residuals is  0.19
## RMSEA index =  0.199  and the 10 % confidence intervals are  0.19 0.208
## BIC =  629.1
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 170  and the fit is  7.91 
## The number of observations was  203  with Chi Square =  1532.35  with prob <  0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012
## The root mean square of the residuals is  0.18 
## The df corrected root mean square of the residuals is  0.19 
## 
## RMSEA index =  0.199  and the 10 % confidence intervals are  0.19 0.208
## BIC =  629.1 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.95   0
## Multiple R square of scores with factors      0.90   0
## Minimum correlation of factor score estimates 0.79  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.88 0.88
## Omega general for total scores and subscales  0.88 0.88
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T3_NOBAGS.1_1r:T3_NOBAGS.16_1r) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.89 
## G.6:                   0.94 
## Omega Hierarchical:    0.89 
## Omega H asymptotic:    0.99 
## Omega Total            0.89 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                     g  F1*   h2   u2 p2
## T3_NOBAGS.1_1r   0.63      0.40 0.60  1
## T3_NOBAGS.1_2r   0.44      0.20 0.80  1
## T3_NOBAGS.2_1r   0.75      0.56 0.44  1
## T3_NOBAGS.2_2r   0.69      0.48 0.52  1
## T3_NOBAGS.3_1r   0.64      0.41 0.59  1
## T3_NOBAGS.3_2r   0.55      0.30 0.70  1
## T3_NOBAGS.4_1r   0.72      0.52 0.48  1
## T3_NOBAGS.4_2r   0.61      0.37 0.63  1
## T3_NOBAGS.5_1r   0.68      0.46 0.54  1
## T3_NOBAGS.6_1r   0.48      0.23 0.77  1
## T3_NOBAGS.7_1r   0.68      0.46 0.54  1
## T3_NOBAGS.8_1r   0.43      0.19 0.81  1
## T3_NOBAGS.9_1r   0.51      0.26 0.74  1
## T3_NOBAGS.10_1r  0.42      0.18 0.82  1
## T3_NOBAGS.11_1r- 0.35      0.12 0.88  1
## T3_NOBAGS.12_1r  0.42      0.18 0.82  1
## T3_NOBAGS.13_1r  0.46      0.21 0.79  1
## T3_NOBAGS.14_1r  0.54      0.29 0.71  1
## T3_NOBAGS.15_1r  0.47      0.22 0.78  1
## T3_NOBAGS.16_1r  0.34      0.12 0.88  1
## 
## With Sums of squares  of:
##   g F1* 
## 6.1 0.0 
## 
## general/max  14727181914176212   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 170  and the fit is  7.67 
## The number of observations was  203  with Chi Square =  1487.54  with prob <  0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000052
## The root mean square of the residuals is  0.17 
## The df corrected root mean square of the residuals is  0.18
## RMSEA index =  0.195  and the 10 % confidence intervals are  0.187 0.205
## BIC =  584.29
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 170  and the fit is  7.67 
## The number of observations was  203  with Chi Square =  1487.54  with prob <  0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000052
## The root mean square of the residuals is  0.17 
## The df corrected root mean square of the residuals is  0.18 
## 
## RMSEA index =  0.195  and the 10 % confidence intervals are  0.187 0.205
## BIC =  584.29 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.95   0
## Multiple R square of scores with factors      0.91   0
## Minimum correlation of factor score estimates 0.82  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.89 0.89
## Omega general for total scores and subscales  0.89 0.89
## Omega group for total scores and subscales    0.00 0.00

Willingness to Help

data <- data %>% 
  mutate(T1_WHS = rowMeans(pick(T1_WHS_1:T1_WHS_6)),
         T2_WHS = rowMeans(pick(T2_WHS_1:T2_WHS_6)),
         T3_WHS = rowMeans(pick(T3_WHS_1:T3_WHS_6)))

# Get alpha & omega
data %>% 
  select(T1_WHS_1:T1_WHS_6) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.65 
## G.6:                   0.66 
## Omega Hierarchical:    0.64 
## Omega H asymptotic:    0.96 
## Omega Total            0.67 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##             g  F1*   h2   u2 p2
## T1_WHS_1 0.72      0.52 0.48  1
## T1_WHS_2 0.57      0.32 0.68  1
## T1_WHS_3 0.24      0.06 0.94  1
## T1_WHS_4 0.46      0.21 0.79  1
## T1_WHS_5 0.25      0.06 0.94  1
## T1_WHS_6 0.67      0.45 0.55  1
## 
## With Sums of squares  of:
##   g F1* 
## 1.6 0.0 
## 
## general/max  14564732675507748   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 9  and the fit is  0.26 
## The number of observations was  203  with Chi Square =  50.63  with prob <  0.000000082
## The root mean square of the residuals is  0.11 
## The df corrected root mean square of the residuals is  0.15
## RMSEA index =  0.151  and the 10 % confidence intervals are  0.112 0.193
## BIC =  2.82
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9  and the fit is  0.26 
## The number of observations was  203  with Chi Square =  50.63  with prob <  0.000000082
## The root mean square of the residuals is  0.11 
## The df corrected root mean square of the residuals is  0.15 
## 
## RMSEA index =  0.151  and the 10 % confidence intervals are  0.112 0.193
## BIC =  2.82 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.86   0
## Multiple R square of scores with factors      0.73   0
## Minimum correlation of factor score estimates 0.46  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.67 0.64
## Omega general for total scores and subscales  0.64 0.64
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T2_WHS_1:T2_WHS_6) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.67 
## G.6:                   0.67 
## Omega Hierarchical:    0.67 
## Omega H asymptotic:    0.98 
## Omega Total            0.68 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##             g  F1*   h2   u2 p2
## T2_WHS_1 0.68      0.47 0.53  1
## T2_WHS_2 0.60      0.36 0.64  1
## T2_WHS_3 0.38      0.14 0.86  1
## T2_WHS_4 0.41      0.16 0.84  1
## T2_WHS_5 0.34      0.11 0.89  1
## T2_WHS_6 0.62      0.39 0.61  1
## 
## With Sums of squares  of:
##   g F1* 
## 1.6 0.0 
## 
## general/max  58885651346390504   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 9  and the fit is  0.26 
## The number of observations was  203  with Chi Square =  51.27  with prob <  0.000000062
## The root mean square of the residuals is  0.11 
## The df corrected root mean square of the residuals is  0.14
## RMSEA index =  0.152  and the 10 % confidence intervals are  0.113 0.194
## BIC =  3.45
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9  and the fit is  0.26 
## The number of observations was  203  with Chi Square =  51.27  with prob <  0.000000062
## The root mean square of the residuals is  0.11 
## The df corrected root mean square of the residuals is  0.14 
## 
## RMSEA index =  0.152  and the 10 % confidence intervals are  0.113 0.194
## BIC =  3.45 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.85   0
## Multiple R square of scores with factors      0.72   0
## Minimum correlation of factor score estimates 0.43  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.68 0.67
## Omega general for total scores and subscales  0.67 0.67
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T3_WHS_1:T3_WHS_6) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.69 
## G.6:                   0.69 
## Omega Hierarchical:    0.69 
## Omega H asymptotic:    0.99 
## Omega Total            0.69 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##             g  F1*   h2   u2 p2
## T3_WHS_1 0.50      0.25 0.75  1
## T3_WHS_2 0.57      0.33 0.67  1
## T3_WHS_3 0.37      0.14 0.86  1
## T3_WHS_4 0.45      0.20 0.80  1
## T3_WHS_5 0.50      0.25 0.75  1
## T3_WHS_6 0.71      0.50 0.50  1
## 
## With Sums of squares  of:
##   g F1* 
## 1.7 0.0 
## 
## general/max  10006094449321290   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 9  and the fit is  0.33 
## The number of observations was  203  with Chi Square =  65.07  with prob <  0.00000000014
## The root mean square of the residuals is  0.12 
## The df corrected root mean square of the residuals is  0.16
## RMSEA index =  0.175  and the 10 % confidence intervals are  0.137 0.217
## BIC =  17.25
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9  and the fit is  0.33 
## The number of observations was  203  with Chi Square =  65.07  with prob <  0.00000000014
## The root mean square of the residuals is  0.12 
## The df corrected root mean square of the residuals is  0.16 
## 
## RMSEA index =  0.175  and the 10 % confidence intervals are  0.137 0.217
## BIC =  17.25 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.85   0
## Multiple R square of scores with factors      0.72   0
## Minimum correlation of factor score estimates 0.44  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.69 0.69
## Omega general for total scores and subscales  0.69 0.69
## Omega group for total scores and subscales    0.00 0.00

Compassionate Love

data <- data %>% 
  mutate(T1_CLS = rowMeans(pick(T1_CLS_1:T1_CLS_21)),
         T2_CLS = rowMeans(pick(T2_CLS_1:T2_CLS_21)),
         T3_CLS = rowMeans(pick(T3_CLS_1:T3_CLS_21)))

# Get alpha & omega
data %>% 
  select(T1_CLS_1:T1_CLS_21) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.94 
## G.6:                   0.95 
## Omega Hierarchical:    0.94 
## Omega H asymptotic:    1 
## Omega Total            0.94 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##              g  F1*   h2   u2 p2
## T1_CLS_1  0.67      0.44 0.56  1
## T1_CLS_2  0.64      0.41 0.59  1
## T1_CLS_3  0.74      0.55 0.45  1
## T1_CLS_4  0.65      0.42 0.58  1
## T1_CLS_5  0.71      0.50 0.50  1
## T1_CLS_6  0.74      0.55 0.45  1
## T1_CLS_7  0.54      0.29 0.71  1
## T1_CLS_8  0.72      0.51 0.49  1
## T1_CLS_9  0.77      0.59 0.41  1
## T1_CLS_10 0.67      0.45 0.55  1
## T1_CLS_11 0.60      0.35 0.65  1
## T1_CLS_12 0.72      0.52 0.48  1
## T1_CLS_13 0.42      0.18 0.82  1
## T1_CLS_14 0.53      0.28 0.72  1
## T1_CLS_15 0.70      0.49 0.51  1
## T1_CLS_16 0.60      0.36 0.64  1
## T1_CLS_17 0.70      0.49 0.51  1
## T1_CLS_18 0.63      0.40 0.60  1
## T1_CLS_19 0.59      0.35 0.65  1
## T1_CLS_20 0.65      0.43 0.57  1
## T1_CLS_21 0.66      0.44 0.56  1
## 
## With Sums of squares  of:
##   g F1* 
##   9   0 
## 
## general/max  54179614793415584   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 189  and the fit is  1.91 
## The number of observations was  203  with Chi Square =  369.45  with prob <  0.000000000000089
## The root mean square of the residuals is  0.06 
## The df corrected root mean square of the residuals is  0.06
## RMSEA index =  0.068  and the 10 % confidence intervals are  0.058 0.079
## BIC =  -634.75
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 189  and the fit is  1.91 
## The number of observations was  203  with Chi Square =  369.45  with prob <  0.000000000000089
## The root mean square of the residuals is  0.06 
## The df corrected root mean square of the residuals is  0.06 
## 
## RMSEA index =  0.068  and the 10 % confidence intervals are  0.058 0.079
## BIC =  -634.75 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.97   0
## Multiple R square of scores with factors      0.94   0
## Minimum correlation of factor score estimates 0.89  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.94 0.94
## Omega general for total scores and subscales  0.94 0.94
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T2_CLS_1:T2_CLS_21) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.94 
## G.6:                   0.96 
## Omega Hierarchical:    0.95 
## Omega H asymptotic:    1 
## Omega Total            0.95 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##              g  F1*   h2   u2 p2
## T2_CLS_1  0.69      0.48 0.52  1
## T2_CLS_2  0.71      0.50 0.50  1
## T2_CLS_3  0.79      0.62 0.38  1
## T2_CLS_4  0.63      0.40 0.60  1
## T2_CLS_5  0.68      0.46 0.54  1
## T2_CLS_6  0.77      0.59 0.41  1
## T2_CLS_7  0.62      0.38 0.62  1
## T2_CLS_8  0.72      0.51 0.49  1
## T2_CLS_9  0.77      0.60 0.40  1
## T2_CLS_10 0.73      0.53 0.47  1
## T2_CLS_11 0.64      0.41 0.59  1
## T2_CLS_12 0.76      0.57 0.43  1
## T2_CLS_13 0.30      0.09 0.91  1
## T2_CLS_14 0.57      0.32 0.68  1
## T2_CLS_15 0.79      0.62 0.38  1
## T2_CLS_16 0.61      0.37 0.63  1
## T2_CLS_17 0.70      0.49 0.51  1
## T2_CLS_18 0.62      0.39 0.61  1
## T2_CLS_19 0.64      0.40 0.60  1
## T2_CLS_20 0.64      0.42 0.58  1
## T2_CLS_21 0.66      0.43 0.57  1
## 
## With Sums of squares  of:
##   g F1* 
## 9.6 0.0 
## 
## general/max  34494890394986312   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 189  and the fit is  2.68 
## The number of observations was  203  with Chi Square =  517.68  with prob <  0.000000000000000000000000000000022
## The root mean square of the residuals is  0.06 
## The df corrected root mean square of the residuals is  0.07
## RMSEA index =  0.092  and the 10 % confidence intervals are  0.083 0.102
## BIC =  -486.51
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 189  and the fit is  2.68 
## The number of observations was  203  with Chi Square =  517.68  with prob <  0.000000000000000000000000000000022
## The root mean square of the residuals is  0.06 
## The df corrected root mean square of the residuals is  0.07 
## 
## RMSEA index =  0.092  and the 10 % confidence intervals are  0.083 0.102
## BIC =  -486.51 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.98   0
## Multiple R square of scores with factors      0.95   0
## Minimum correlation of factor score estimates 0.90  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.95 0.95
## Omega general for total scores and subscales  0.95 0.95
## Omega group for total scores and subscales    0.00 0.00
data %>% 
  select(T3_CLS_1:T3_CLS_21) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.95 
## G.6:                   0.96 
## Omega Hierarchical:    0.95 
## Omega H asymptotic:    1 
## Omega Total            0.95 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##              g  F1*   h2   u2 p2
## T3_CLS_1  0.73      0.53 0.47  1
## T3_CLS_2  0.76      0.57 0.43  1
## T3_CLS_3  0.77      0.59 0.41  1
## T3_CLS_4  0.71      0.50 0.50  1
## T3_CLS_5  0.69      0.48 0.52  1
## T3_CLS_6  0.79      0.63 0.37  1
## T3_CLS_7  0.65      0.42 0.58  1
## T3_CLS_8  0.71      0.50 0.50  1
## T3_CLS_9  0.77      0.59 0.41  1
## T3_CLS_10 0.73      0.54 0.46  1
## T3_CLS_11 0.66      0.43 0.57  1
## T3_CLS_12 0.79      0.62 0.38  1
## T3_CLS_13 0.46      0.21 0.79  1
## T3_CLS_14 0.55      0.30 0.70  1
## T3_CLS_15 0.80      0.64 0.36  1
## T3_CLS_16 0.61      0.38 0.62  1
## T3_CLS_17 0.74      0.55 0.45  1
## T3_CLS_18 0.69      0.48 0.52  1
## T3_CLS_19 0.61      0.37 0.63  1
## T3_CLS_20 0.63      0.39 0.61  1
## T3_CLS_21 0.68      0.46 0.54  1
## 
## With Sums of squares  of:
##   g F1* 
##  10   0 
## 
## general/max  36644823991854128   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 189  and the fit is  2.78 
## The number of observations was  203  with Chi Square =  537.38  with prob <  0.000000000000000000000000000000000038
## The root mean square of the residuals is  0.06 
## The df corrected root mean square of the residuals is  0.07
## RMSEA index =  0.095  and the 10 % confidence intervals are  0.086 0.105
## BIC =  -466.81
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 189  and the fit is  2.78 
## The number of observations was  203  with Chi Square =  537.38  with prob <  0.000000000000000000000000000000000038
## The root mean square of the residuals is  0.06 
## The df corrected root mean square of the residuals is  0.07 
## 
## RMSEA index =  0.095  and the 10 % confidence intervals are  0.086 0.105
## BIC =  -466.81 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.98   0
## Multiple R square of scores with factors      0.96   0
## Minimum correlation of factor score estimates 0.91  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.95 0.95
## Omega general for total scores and subscales  0.95 0.95
## Omega group for total scores and subscales    0.00 0.00

Situational Self-Control

data <- data %>% 
  mutate(across(contains("SMS5"), .names = "{col}r"))

# Reverse code SMS5 (items 3 et 5)
data <- data %>%
  mutate(across(ends_with(paste0("SMS5_", c(3, 5))), 
                ~nice_reverse(.x, 5), .names = "{col}r"))

# Get mean SMS5
data <- data %>% 
  mutate(T2_SMS5 = rowMeans(pick(T2_SMS5_1r:T2_SMS5_6r)))

# Get alpha & omega
data %>% 
  select(T2_SMS5_1r:T2_SMS5_6r) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.65 
## G.6:                   0.65 
## Omega Hierarchical:    0.67 
## Omega H asymptotic:    1 
## Omega Total            0.67 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##               g  F1*   h2   u2 p2
## T2_SMS5_1r 0.63      0.40 0.60  1
## T2_SMS5_2r           0.00 1.00  1
## T2_SMS5_3r 0.50      0.25 0.75  1
## T2_SMS5_4r 0.58      0.34 0.66  1
## T2_SMS5_5r 0.59      0.35 0.65  1
## T2_SMS5_6r 0.62      0.39 0.61  1
## 
## With Sums of squares  of:
##   g F1* 
## 1.7 0.0 
## 
## general/max  11338623069083538   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 9  and the fit is  0.15 
## The number of observations was  203  with Chi Square =  30.12  with prob <  0.00042
## The root mean square of the residuals is  0.07 
## The df corrected root mean square of the residuals is  0.09
## RMSEA index =  0.107  and the 10 % confidence intervals are  0.067 0.151
## BIC =  -17.7
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9  and the fit is  0.15 
## The number of observations was  203  with Chi Square =  30.12  with prob <  0.00042
## The root mean square of the residuals is  0.07 
## The df corrected root mean square of the residuals is  0.09 
## 
## RMSEA index =  0.107  and the 10 % confidence intervals are  0.067 0.151
## BIC =  -17.7 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.85   0
## Multiple R square of scores with factors      0.73   0
## Minimum correlation of factor score estimates 0.46  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.67 0.67
## Omega general for total scores and subscales  0.67 0.67
## Omega group for total scores and subscales    0.00 0.00

Positive and Negative Affect

# Get mean of PANAS
# Positive affect = 1, 3, 5, 7, 9
# Negative affect = 2, 4, 6, 8, 10
data <- data %>% mutate(
  T2_PANAS_pos = rowMeans(pick(paste0("T2_PANAS_", seq(1, 9, 2)))),
  T2_PANAS_neg = rowMeans(pick(paste0("T2_PANAS_", seq(2, 10, 2)))))

# Get alpha & omega
# Positive affect
data %>% 
  select(all_of(paste0("T2_PANAS_", seq(1, 9, 2)))) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in cov2cor(t(w) %*% r %*% w): diag(.) had 0 or NA entries; non-finite
## result is doubtful
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.86 
## G.6:                   0.84 
## Omega Hierarchical:    0.86 
## Omega H asymptotic:    1 
## Omega Total            0.86 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##               g  F1*   h2   u2 p2
## T2_PANAS_1 0.73      0.53 0.47  1
## T2_PANAS_3 0.70      0.49 0.51  1
## T2_PANAS_5 0.79      0.63 0.37  1
## T2_PANAS_7 0.67      0.45 0.55  1
## T2_PANAS_9 0.84      0.70 0.30  1
## 
## With Sums of squares  of:
##   g F1* 
## 2.8 0.0 
## 
## general/max  Inf   max/min =   NaN
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 5  and the fit is  0.05 
## The number of observations was  203  with Chi Square =  10.65  with prob <  0.059
## The root mean square of the residuals is  0.03 
## The df corrected root mean square of the residuals is  0.04
## RMSEA index =  0.074  and the 10 % confidence intervals are  0 0.138
## BIC =  -15.91
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 5  and the fit is  0.05 
## The number of observations was  203  with Chi Square =  10.65  with prob <  0.059
## The root mean square of the residuals is  0.03 
## The df corrected root mean square of the residuals is  0.04 
## 
## RMSEA index =  0.074  and the 10 % confidence intervals are  0 0.138
## BIC =  -15.91 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.94   0
## Multiple R square of scores with factors      0.88   0
## Minimum correlation of factor score estimates 0.75  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.86 0.86
## Omega general for total scores and subscales  0.86 0.86
## Omega group for total scores and subscales    0.00 0.00
# Negative affect
data %>% 
  select(all_of(paste0("T2_PANAS_", seq(2, 10, 2)))) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor

## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## diag(.) had 0 or NA entries; non-finite result is doubtful
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.8 
## G.6:                   0.8 
## Omega Hierarchical:    0.82 
## Omega H asymptotic:    1 
## Omega Total            0.82 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                g  F1*   h2   u2 p2
## T2_PANAS_2  0.85      0.73 0.27  1
## T2_PANAS_4  0.34      0.11 0.89  1
## T2_PANAS_6  0.79      0.62 0.38  1
## T2_PANAS_8  0.67      0.45 0.55  1
## T2_PANAS_10 0.73      0.54 0.46  1
## 
## With Sums of squares  of:
##   g F1* 
## 2.5 0.0 
## 
## general/max  Inf   max/min =   NaN
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 5  and the fit is  0.16 
## The number of observations was  203  with Chi Square =  31.84  with prob <  0.0000064
## The root mean square of the residuals is  0.07 
## The df corrected root mean square of the residuals is  0.1
## RMSEA index =  0.163  and the 10 % confidence intervals are  0.112 0.219
## BIC =  5.28
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 5  and the fit is  0.16 
## The number of observations was  203  with Chi Square =  31.84  with prob <  0.0000064
## The root mean square of the residuals is  0.07 
## The df corrected root mean square of the residuals is  0.1 
## 
## RMSEA index =  0.163  and the 10 % confidence intervals are  0.112 0.219
## BIC =  5.28 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.93   0
## Multiple R square of scores with factors      0.87   0
## Minimum correlation of factor score estimates 0.73  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.82 0.82
## Omega general for total scores and subscales  0.82 0.82
## Omega group for total scores and subscales    0.00 0.00

Charity

data <- data %>% mutate(
  T2_Charity = rowMeans(pick(contains("charity") & ends_with("1_1"))),
  T2_Familiarity = rowMeans(pick(contains("charity") & ends_with("2_1"))))

# Get alpha & omega
data %>% 
  select(contains("charity") & ends_with("1_1")) %>% 
  omega(nfactors = 1)
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in cov2cor(t(w) %*% r %*% w): diag(.) had 0 or NA entries; non-finite
## result is doubtful
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.98 
## G.6:                   0.99 
## Omega Hierarchical:    0.98 
## Omega H asymptotic:    1 
## Omega Total            0.98 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##                            g  F1*   h2   u2 p2
## T2_charity.moisson1_1   0.84      0.71 0.29  1
## T2_charity.bonneau1_1   0.85      0.72 0.28  1
## T2_charity.AMDI1_1      0.86      0.74 0.26  1
## T2_charity.CAAM1_1      0.86      0.74 0.26  1
## T2_charity.centrefem1_1 0.80      0.64 0.36  1
## T2_charity.LGBTQ1_1     0.74      0.55 0.45  1
## T2_charity.equiterre1_1 0.87      0.75 0.25  1
## T2_charity.dejeuner1_1  0.80      0.65 0.35  1
## T2_charity.cancer1_1    0.87      0.75 0.25  1
## T2_charity.parkinson1_1 0.87      0.76 0.24  1
## T2_charity.OXFAM1_1     0.81      0.66 0.34  1
## T2_charity.papillon1_1  0.85      0.72 0.28  1
## T2_charity.croix1_1     0.86      0.74 0.26  1
## T2_charity.centraide1_1 0.87      0.76 0.24  1
## T2_charity.autisme1_1   0.82      0.68 0.32  1
## T2_charity.suzuki1_1    0.78      0.62 0.38  1
## T2_charity.conserv1_1   0.81      0.65 0.35  1
## T2_charity.coeur1_1     0.86      0.73 0.27  1
## T2_charity.UNICEF1_1    0.82      0.68 0.32  1
## T2_charity.amnistie1_1  0.81      0.66 0.34  1
## T2_charity.green1_1     0.77      0.59 0.41  1
## T2_charity.WWF1_1       0.79      0.62 0.38  1
## T2_charity.MSF1_1       0.86      0.74 0.26  1
## T2_charity.armee1_1     0.75      0.57 0.43  1
## 
## With Sums of squares  of:
##   g F1* 
##  16   0 
## 
## general/max  Inf   max/min =   NaN
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 252  and the fit is  4.94 
## The number of observations was  203  with Chi Square =  950.72  with prob <  0.0000000000000000000000000000000000000000000000000000000000000000000000000000000011
## The root mean square of the residuals is  0.05 
## The df corrected root mean square of the residuals is  0.05
## RMSEA index =  0.117  and the 10 % confidence intervals are  0.109 0.125
## BIC =  -388.21
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 252  and the fit is  4.94 
## The number of observations was  203  with Chi Square =  950.72  with prob <  0.0000000000000000000000000000000000000000000000000000000000000000000000000000000011
## The root mean square of the residuals is  0.05 
## The df corrected root mean square of the residuals is  0.05 
## 
## RMSEA index =  0.117  and the 10 % confidence intervals are  0.109 0.125
## BIC =  -388.21 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.99   0
## Multiple R square of scores with factors      0.98   0
## Minimum correlation of factor score estimates 0.96  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.98 0.98
## Omega general for total scores and subscales  0.98 0.98
## Omega group for total scores and subscales    0.00 0.00

Intensity * Duration

# Create new variable blastintensity.duration
data <- data %>% 
  mutate(T1_blastintensity.duration = T1_blastintensity * T1_blastduration,
         T2_blastintensity.duration = T2_blastintensity * T2_blastduration,
         T3_blastintensity.duration = T3_blastintensity * T3_blastduration)

Assumptions

In this section, we: (a) test assumptions of normality, (b) transform variables violating assumptions, (c) test assumptions of homoscedasticity, and (d) identify and winsorize outliers.

…

At this stage, we define a list of our relevant variables.

# Make list of DVs
col.list <- c("T1_blastintensity", "T1_blastduration", "T1_blastintensity.duration",
              "T2_blastintensity", "T2_blastduration", "T2_blastintensity.duration",
              "T3_blastintensity", "T3_blastduration", "T3_blastintensity.duration",
              "T1_BSCS", "T1_BAQ", "T1_attitude", "T2_attitude", "T3_attitude",
              "T1_dehumanization", "T2_dehumanization", "T3_dehumanization", 
              "T1_NOBAGS", "T1_WHS", "T2_WHS", "T3_WHS", "T1_CLS", "T2_CLS", "T3_CLS", 
              "T2_SMS5", "T2_PANAS_pos", "T2_PANAS_neg", "T2_Charity", 
              "T2_Familiarity", "T2_memory.altruistic", "T2_memory.aggressive")

Normality

lapply(col.list, function(x) 
  nice_normality(data, 
                 variable = x, 
                 title = x,
                 group = "T1_Group",
                 shapiro = TRUE,
                 histogram = TRUE))
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

Several variables are clearly skewed. Let’s apply transformations. But first, let’s deal with the working memory task, SOPT (Self-Ordered Pointing Task). It is clearly problematic.

Transformation

The function below transforms variables according to the best possible transformation (via the bestNormalize package), and also standardizes the variables.

Note. The I(x) transformations above are actually not transformations, but a shorthand function for passing the data “as is”. Suggesting the package estimated the various attempted transformations did not improve normality in those cases, so no transformation is used. This only appears when standardize is set to FALSE. When set to TRUE, for those variables, it is actually center_scale(x), suggesting that the data are only CENTERED because they need no transformation (no need to be scaled), only to be centered.

Let’s check if normality was corrected.

Looks rather reasonable now, though not perfect (fortunately t-tests are quite robust against violations of normality).

We can now resume with the next step: checking variance.

Homoscedasticity

# Plotting variance
plots(lapply(col.list, function(x) {
  nice_varplot(data, x, group = "T1_Group")
  }),
  n_columns = 2)

Variance looks good. No group has four times the variance of any other group. We can now resume with checking outliers.

Outliers

We check outliers visually with the plot_outliers function, which draws red lines at +/- 3 median absolute deviations.

plots(lapply(col.list, function(x) {
  plot_outliers(data, x, group = "T1_Group", ytitle = x, binwidth = 0.3)
  }),
  n_columns = 2)

There are some outliers, but nothing unreasonable. Let’s still check with the 3 median absolute deviations (MAD) method.

data %>% 
  as.data.frame %>% 
  filter(T1_Group == "Waitlist") %>% 
  find_mad(col.list, criteria = 3)
## 39 outlier(s) based on 3 median absolute deviations for variable(s): 
##  T1_blastintensity, T1_blastduration, T1_blastintensity.duration, T2_blastintensity, T2_blastduration, T2_blastintensity.duration, T3_blastintensity, T3_blastduration, T3_blastintensity.duration, T1_BSCS, T1_BAQ, T1_attitude, T2_attitude, T3_attitude, T1_dehumanization, T2_dehumanization, T3_dehumanization, T1_NOBAGS, T1_WHS, T2_WHS, T3_WHS, T1_CLS, T2_CLS, T3_CLS, T2_SMS5, T2_PANAS_pos, T2_PANAS_neg, T2_Charity, T2_Familiarity, T2_memory.altruistic, T2_memory.aggressive 
## 
## The following participants were considered outliers for more than one variable: 
## 
##    Row n
## 1    6 6
## 2    7 5
## 3   16 2
## 4   20 3
## 5   26 2
## 6   37 3
## 7   43 2
## 8   49 2
## 9   50 2
## 10  51 3
## 11  66 2
## 12  69 3
## 13  70 2
## 14  83 2
## 15  87 4
## 
## Outliers per variable: 
## 
## $T1_blastintensity.duration
##   Row T1_blastintensity.duration_mad
## 1  37                       4.279092
## 
## $T3_blastintensity.duration
##   Row T3_blastintensity.duration_mad
## 1   4                       3.432722
## 2   9                       3.353952
## 
## $T1_attitude
##   Row T1_attitude_mad
## 1   6       -5.178348
## 2  15       -3.054790
## 3  26       -3.220149
## 4  50       -3.098306
## 
## $T2_attitude
##   Row T2_attitude_mad
## 1   6       -5.750247
## 2  26       -3.052284
## 3  66       -3.052284
## 
## $T3_attitude
##   Row T3_attitude_mad
## 1   6       -6.397443
## 
## $T1_dehumanization
##    Row T1_dehumanization_mad
## 1    5             -8.408651
## 2    6            -13.040155
## 3    7             -5.395926
## 4   16             -7.329466
## 5   20             -6.542560
## 6   41             -3.934529
## 7   43             -6.902289
## 8   50             -3.709699
## 9   51             -5.440892
## 10  58             -5.845587
## 11  64             -3.507352
## 12  69             -6.025451
## 13  79             -3.147624
## 14  83             -4.721435
## 15  87             -5.575790
## 16  88             -4.609020
## 
## $T2_dehumanization
##   Row T2_dehumanization_mad
## 1   6            -13.799716
## 2   7             -3.919338
## 3  20             -5.541762
## 4  37             -3.208388
## 5  51             -3.500060
## 6  69             -3.554749
## 7  70             -4.265698
## 8  87             -5.049566
## 
## $T3_dehumanization
##    Row T3_dehumanization_mad
## 1    6            -13.267943
## 2    7             -3.594326
## 3   16             -4.393065
## 4   20             -7.304025
## 5   37             -3.008584
## 6   43             -3.416828
## 7   51             -5.830795
## 8   69             -5.475800
## 9   70             -4.144568
## 10  83             -3.328079
## 
## $T3_WHS
##   Row T3_WHS_mad
## 1  22  -3.083386
## 
## $T2_CLS
##   Row T2_CLS_mad
## 1  87  -3.480372
## 
## $T3_CLS
##   Row T3_CLS_mad
## 1  10   -3.06115
## 
## $T2_PANAS_neg
##   Row T2_PANAS_neg_mad
## 1   7         3.035208
## 2  44         5.058681
## 
## $T2_memory.altruistic
##   Row T2_memory.altruistic_mad
## 1   8                 5.232062
## 2  18                 3.080544
## 3  49                 3.726311
## 4  56                 8.054610
## 5  66                 9.220174
## 6  77                 3.873390
## 7  87                 3.803484
## 8  90                 4.521580
## 
## $T2_memory.aggressive
##   Row T2_memory.aggressive_mad
## 1   1                 3.045000
## 2   2                 5.060588
## 3   7                 5.035664
## 4  30                 4.291231
## 5  34                 5.927863
## 6  48                 3.351726
## 7  49                 3.652856
## 8  68                 4.764036
## 9  84                 4.659505
data %>% 
  as.data.frame %>% 
  filter(T1_Group == "Reflection") %>% 
  find_mad(col.list, criteria = 3)
## 33 outlier(s) based on 3 median absolute deviations for variable(s): 
##  T1_blastintensity, T1_blastduration, T1_blastintensity.duration, T2_blastintensity, T2_blastduration, T2_blastintensity.duration, T3_blastintensity, T3_blastduration, T3_blastintensity.duration, T1_BSCS, T1_BAQ, T1_attitude, T2_attitude, T3_attitude, T1_dehumanization, T2_dehumanization, T3_dehumanization, T1_NOBAGS, T1_WHS, T2_WHS, T3_WHS, T1_CLS, T2_CLS, T3_CLS, T2_SMS5, T2_PANAS_pos, T2_PANAS_neg, T2_Charity, T2_Familiarity, T2_memory.altruistic, T2_memory.aggressive 
## 
## The following participants were considered outliers for more than one variable: 
## 
##    Row n
## 1    2 2
## 2    4 2
## 3    5 3
## 4    8 3
## 5   22 6
## 6   25 6
## 7   29 4
## 8   31 2
## 9   36 2
## 10  37 5
## 11  39 5
## 12  43 4
## 13  50 2
## 14  53 5
## 
## Outliers per variable: 
## 
## $T1_blastduration
##   Row T1_blastduration_mad
## 1  50             -3.01203
## 
## $T1_blastintensity.duration
##   Row T1_blastintensity.duration_mad
## 1   8                       4.013772
## 2  41                       3.607649
## 
## $T2_blastintensity.duration
##   Row T2_blastintensity.duration_mad
## 1   5                       3.228501
## 2   8                       7.975666
## 3  43                       3.268278
## 
## $T3_blastintensity.duration
##   Row T3_blastintensity.duration_mad
## 1   5                       3.777493
## 2   8                       7.395048
## 3  36                       3.236847
## 4  37                       3.781831
## 5  42                       5.320105
## 
## $T2_attitude
##    Row T2_attitude_mad
## 1    2       -3.489757
## 2   18      -11.524994
## 3   22       -5.190646
## 4   23       -3.357791
## 5   25       -5.410589
## 6   29       -6.407662
## 7   37       -5.615869
## 8   39       -3.123185
## 9   43       -3.357791
## 10  53       -4.105596
## 
## $T3_attitude
##   Row T3_attitude_mad
## 1  11       -7.503710
## 2  22       -6.070417
## 3  25       -5.817483
## 4  29       -6.306489
## 5  37       -5.261028
## 6  39       -5.328477
## 7  43       -4.974369
## 8  48       -3.878322
## 9  53       -7.874680
## 
## $T1_dehumanization
##   Row T1_dehumanization_mad
## 1  25            -10.791852
## 2  34             -3.736679
## 3  39             -4.910293
## 4  43             -5.935519
## 
## $T2_dehumanization
##   Row T2_dehumanization_mad
## 1  22             -4.903029
## 2  25             -4.903029
## 3  29             -4.928971
## 4  37             -4.474987
## 5  38             -4.124770
## 6  39             -4.215567
## 7  53             -3.891293
## 
## $T3_dehumanization
##   Row T3_dehumanization_mad
## 1  22             -3.839409
## 2  25             -3.320570
## 3  29             -3.663004
## 4  39             -3.901670
## 5  53             -3.548859
## 
## $T1_WHS
##   Row T1_WHS_mad
## 1  52  -3.597284
## 
## $T2_SMS5
##   Row T2_SMS5_mad
## 1  22    3.035208
## 2  50    3.709699
## 
## $T2_PANAS_neg
##   Row T2_PANAS_neg_mad
## 1   2         3.372454
## 2  12         3.372454
## 3  24         3.372454
## 4  53         4.046945
## 
## $T2_Familiarity
##   Row T2_Familiarity_mad
## 1  20           3.065867
## 2  36           3.065867
## 3  37           3.556406
## 
## $T2_memory.altruistic
##   Row T2_memory.altruistic_mad
## 1   3                 3.003502
## 2   4                 3.129885
## 3  30                 3.723668
## 4  31                 3.455823
## 5  49                 3.000397
## 
## $T2_memory.aggressive
##   Row T2_memory.aggressive_mad
## 1   1                 4.158773
## 2   4                 6.946257
## 3   5                 6.538862
## 4  15                 3.383314
## 5  16                 3.669782
## 6  22                11.217445
## 7  25                 3.247418
## 8  31                 5.680633
## 9  44                10.604004
data %>% 
  as.data.frame %>% 
  filter(T1_Group == "Meditation") %>% 
  find_mad(col.list, criteria = 3)
## 36 outlier(s) based on 3 median absolute deviations for variable(s): 
##  T1_blastintensity, T1_blastduration, T1_blastintensity.duration, T2_blastintensity, T2_blastduration, T2_blastintensity.duration, T3_blastintensity, T3_blastduration, T3_blastintensity.duration, T1_BSCS, T1_BAQ, T1_attitude, T2_attitude, T3_attitude, T1_dehumanization, T2_dehumanization, T3_dehumanization, T1_NOBAGS, T1_WHS, T2_WHS, T3_WHS, T1_CLS, T2_CLS, T3_CLS, T2_SMS5, T2_PANAS_pos, T2_PANAS_neg, T2_Charity, T2_Familiarity, T2_memory.altruistic, T2_memory.aggressive 
## 
## The following participants were considered outliers for more than one variable: 
## 
##    Row n
## 1    1 2
## 2    3 4
## 3    5 3
## 4   10 2
## 5   11 3
## 6   16 2
## 7   20 3
## 8   26 2
## 9   35 4
## 10  38 4
## 11  41 3
## 12  42 2
## 13  48 5
## 14  50 3
## 15  53 5
## 16  55 4
## 17  56 2
## 18  58 4
## 
## Outliers per variable: 
## 
## $T1_blastintensity
##   Row T1_blastintensity_mad
## 1  58              3.102657
## 
## $T1_blastduration
##   Row T1_blastduration_mad
## 1  19            -3.157726
## 2  27            -3.157726
## 3  38            -3.157726
## 4  58             3.268636
## 
## $T1_blastintensity.duration
##   Row T1_blastintensity.duration_mad
## 1  37                       3.001374
## 2  53                       3.041189
## 3  58                       5.630293
## 
## $T2_blastintensity.duration
##   Row T2_blastintensity.duration_mad
## 1   6                       3.004049
## 
## $T1_attitude
##   Row T1_attitude_mad
## 1   3       -6.936605
## 2   5       -5.715422
## 3  11       -4.437439
## 4  35       -5.218429
## 5  38       -3.869447
## 6  41       -5.289428
## 7  48       -5.630223
## 8  55       -4.153443
## 
## $T2_attitude
##    Row T2_attitude_mad
## 1    1       -3.019149
## 2    5       -6.552196
## 3   20       -6.552196
## 4   35       -6.231010
## 5   41       -5.733171
## 6   42       -9.619523
## 7   48       -6.423722
## 8   49       -3.982707
## 9   50       -4.769613
## 10  55       -5.990120
## 
## $T3_attitude
##    Row T3_attitude_mad
## 1    1       -3.495088
## 2    3      -10.460739
## 3    5      -10.362631
## 4   16       -6.610009
## 5   35      -10.362631
## 6   41       -5.825147
## 7   45       -3.568669
## 8   48       -9.234392
## 9   50       -4.083735
## 10  53       -6.094944
## 11  55      -10.411685
## 12  56       -3.372454
## 13  58      -10.779589
## 
## $T1_dehumanization
##   Row T1_dehumanization_mad
## 1  16             -3.912046
## 2  31             -8.129389
## 3  48             -6.027817
## 4  53             -6.084617
## 5  57             -3.528652
## 
## $T2_dehumanization
##   Row T2_dehumanization_mad
## 1  42             -9.011510
## 2  50             -3.443040
## 3  53             -3.097952
## 
## $T3_dehumanization
##   Row T3_dehumanization_mad
## 1   3             -5.039945
## 2   8             -5.658228
## 3  35             -7.381927
## 4  53             -5.377190
## 
## $T1_CLS
##   Row T1_CLS_mad
## 1  38  -3.460431
## 2  48  -4.398853
## 
## $T3_CLS
##   Row T3_CLS_mad
## 1  38  -3.059297
## 
## $T2_SMS5
##   Row T2_SMS5_mad
## 1  21    3.102657
## 
## $T2_PANAS_neg
##   Row T2_PANAS_neg_mad
## 1  13         3.709699
## 2  40         5.058681
## 3  54         3.372454
## 4  56         3.709699
## 
## $T2_memory.altruistic
##   Row T2_memory.altruistic_mad
## 1   3                 3.917226
## 2  10                 3.177224
## 3  11                 7.148862
## 4  20                 5.338032
## 5  23                 3.700884
## 6  26                 4.172745
## 7  28                 3.457989
## 8  55                 9.298788
## 
## $T2_memory.aggressive
##   Row T2_memory.aggressive_mad
## 1   2                 7.108767
## 2  10                20.601909
## 3  11                 3.306058
## 4  18                22.846607
## 5  20                14.980875
## 6  26                 3.071248
## 7  34                 3.698610

After our transformations, there are 12 outliers in the waitlist group, 3 in the reflection group, and 12 in the meditation group.

Multivariate outliers

For multivariate outliers, it is recommended to use the Minimum Covariance Determinant, a robust version of the Mahalanobis distance (MCD, Leys et al., 2019).

Leys, C., Delacre, M., Mora, Y. L., Lakens, D., & Ley, C. (2019). How to classify, detect, and manage univariate and multivariate outliers, with emphasis on pre-registration. International Review of Social Psychology, 32(1).

data.na <- na.omit(data[col.list])
x <- check_outliers(data.na[-c(15:16)], method = "mcd", threshold = 500)
# We have to exclude two variables that are too large following the transformations
# Otherwise we get an error:
# Error in solve.default(cov, ...) : 
#  system is computationally singular: reciprocal condition number = 8.99059e-20
x
## OK: No outliers detected.
## - Based on the following method and threshold: mcd (500).
## - For variables: T1_blastintensity, T1_blastduration, T1_blastintensity.duration, T2_blastintensity, T2_blastduration, T2_blastintensity.duration, T3_blastintensity, T3_blastduration, T3_blastintensity.duration, T1_BSCS, T1_BAQ, T1_attitude, T2_attitude, T3_attitude, T3_dehumanization, T1_NOBAGS, T1_WHS, T2_WHS, T3_WHS, T1_CLS, T2_CLS, T3_CLS, T2_SMS5, T2_PANAS_pos, T2_PANAS_neg, T2_Charity, T2_Familiarity, T2_memory.altruistic, T2_memory.aggressive
plot(x)

There are 4 multivariate outliers according to the MCD method using an artificially high threshold. However, we did not mention in the preregistration that we would exclude multivariate outliers, so we will not at this time.

Winsorization

Visual assessment and the MAD method confirm we have some outlier values. We could ignore them but because they could have disproportionate influence on the models, one recommendation is to winsorize them by bringing the values at 3 SD. Instead of using the standard deviation around the mean, however, we use the absolute deviation around the median, as it is more robust to extreme observations. For a discussion, see:

Leys, C., Klein, O., Bernard, P., & Licata, L. (2013). Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median. Journal of Experimental Social Psychology, 49(4), 764–766. https://doi.org/10.1016/j.jesp.2013.03.013

Outliers are still present but were brought back within reasonable limits, where applicable.

Standardization

We can now standardize our variables.

We are now ready to move to the analyses proper.

# If we decide to only center variables instead of standardizing them
# (as in the preregistration), then let's do this instead
data <- data %>%
  mutate(across(all_of(gsub(".t.w.s", "", col.list)), 
                \(x) standardize(x, center = TRUE, scale = FALSE), 
                .names = "{col}"))

Contrasts (Group Differences)

In this section, we perform the planned contrasts.

…

Time 2

# Specify the order of factor levels for "Group". 
# Otherwise R will alphabetize them.
data$T1_Group <- factor(data$T1_Group, levels = c("Meditation", "Reflection", "Waitlist"))

# Define our dependent variables
DV <- data %>% select(T2_NOBAGS:T2_Charity) %>% names

# First column (which variable)
Variable <- rep(DV, each = 3)

# Second column (which comparison)
Comparison <- rep(c("MeditationvsCTR", 
                    "ReflectionvsCTR", 
                    "MeditationvsReflection"), 
                  length(DV))
# 14 == number of DV

# Make list of all formulas
formulas <- c(
  "T2_NOBAGS ~ T1_Group * T1_NOBAGS",
  "T2_attitude ~ T1_Group * T1_attitude",
  "T2_dehumanization ~ T1_Group * T1_dehumanization",
  "T2_IAT ~ T1_Group * T1_IAT",
  "T2_SMS5 ~ T1_Group",
  "T2_PANAS_pos ~ T1_Group",
  "T2_PANAS_neg ~ T1_Group",
  "T2_blastintensity ~ T1_Group * T1_blastintensity",
  "T2_blastduration ~ T1_Group * T1_blastduration",
  "T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration",
  "T2_memory.altruistic ~ T1_Group",
  "T2_memory.aggressive ~ T1_Group",
  "T2_WHS ~ T1_Group * T1_WHS",
  "T2_CLS ~ T1_Group * T1_CLS",
  "T2_Charity ~ T1_Group * T2_Familiarity" 
)

# Make list of all models
models.list <- sapply(formulas, lm, data = data, simplify = FALSE, USE.NAMES = TRUE)

## Attempt with nice_lm_contrasts 
x <- nice_lm_contrasts(models.list, group = "T1_Group", data = data)

x.table <- nice_table(x, highlight = TRUE)
x.table

Dependent Variable

Comparison

df

t

p

d

95% CI

T2_NOBAGS

Meditation - Waitlist

197

-1.47

.144

-0.30

[-0.64, 0.03]

Reflection - Waitlist

197

-1.40

.162

-0.21

[-0.59, 0.13]

Meditation - Reflection

197

-0.03

.979

-0.09

[-0.50, 0.29]

T2_attitude

Meditation - Waitlist

197

0.82

.411

0.17

[-0.18, 0.51]

Reflection - Waitlist

197

1.25

.213

0.08

[-0.28, 0.43]

Meditation - Reflection

197

-0.41

.684

0.09

[-0.29, 0.50]

T2_dehumanization

Meditation - Waitlist

195

0.82

.415

0.13

[-0.20, 0.39]

Reflection - Waitlist

195

-0.65

.518

-0.15

[-0.53, 0.22]

Meditation - Reflection

195

1.31

.191

0.28

[-0.10, 0.65]

T2_IAT

Meditation - Waitlist

197

0.51

.611

0.08

[-0.23, 0.39]

Reflection - Waitlist

197

1.14

.256

0.38

[0.01, 0.75]

Meditation - Reflection

197

-0.59

.554

-0.31

[-0.68, 0.06]

T2_SMS5

Meditation - Waitlist

200

-1.50

.136

-0.25

[-0.55, 0.08]

Reflection - Waitlist

200

-0.16

.876

-0.03

[-0.38, 0.32]

Meditation - Reflection

200

-1.18

.239

-0.22

[-0.59, 0.15]

T2_PANAS_pos

Meditation - Waitlist

200

2.76

.006**

0.46

[0.13, 0.82]

Reflection - Waitlist

200

1.17

.244

0.20

[-0.15, 0.53]

Meditation - Reflection

200

1.38

.170

0.26

[-0.12, 0.65]

T2_PANAS_neg

Meditation - Waitlist

200

1.60

.111

0.27

[-0.07, 0.62]

Reflection - Waitlist

200

0.44

.662

0.08

[-0.27, 0.40]

Meditation - Reflection

200

1.01

.312

0.19

[-0.19, 0.57]

T2_blastintensity

Meditation - Waitlist

197

-1.35

.179

-0.03

[-0.37, 0.29]

Reflection - Waitlist

197

-4.16

< .001***

-0.19

[-0.57, 0.19]

Meditation - Reflection

197

2.62

.009**

0.16

[-0.22, 0.57]

T2_blastduration

Meditation - Waitlist

197

-0.83

.408

-0.02

[-0.34, 0.30]

Reflection - Waitlist

197

-4.29

< .001***

-0.33

[-0.68, 0.05]

Meditation - Reflection

197

3.18

.002**

0.31

[-0.08, 0.70]

T2_blastintensity.duration

Meditation - Waitlist

197

-1.06

.291

-0.04

[-0.36, 0.27]

Reflection - Waitlist

197

-3.20

.002**

-0.14

[-0.52, 0.25]

Meditation - Reflection

197

1.99

.048*

0.10

[-0.30, 0.51]

T2_memory.altruistic

Meditation - Waitlist

200

-0.52

.603

-0.09

[-0.42, 0.30]

Reflection - Waitlist

200

-2.44

.015*

-0.42

[-0.65, -0.15]

Meditation - Reflection

200

1.76

.080

0.33

[0.06, 0.68]

T2_memory.aggressive

Meditation - Waitlist

200

1.36

.176

0.23

[-0.16, 0.59]

Reflection - Waitlist

200

0.36

.717

0.06

[-0.21, 0.36]

Meditation - Reflection

200

0.87

.386

0.17

[-0.27, 0.56]

T2_WHS

Meditation - Waitlist

197

0.64

.524

0.07

[-0.25, 0.42]

Reflection - Waitlist

197

1.24

.218

-0.01

[-0.35, 0.31]

Meditation - Reflection

197

-0.57

.573

0.08

[-0.31, 0.45]

T2_CLS

Meditation - Waitlist

197

2.56

.011*

0.40

[0.04, 0.72]

Reflection - Waitlist

197

3.46

.001***

0.43

[0.08, 0.77]

Meditation - Reflection

197

-0.88

.382

-0.03

[-0.43, 0.37]

T2_Charity

Meditation - Waitlist

189

-1.55

.124

-0.19

[-0.54, 0.15]

Reflection - Waitlist

189

-0.93

.355

-0.09

[-0.43, 0.25]

Meditation - Reflection

189

-0.52

.601

-0.10

[-0.49, 0.34]

Time 3

# Make list of all formulas
formulas2 <- c(
  "T3_NOBAGS ~ T1_Group * T1_NOBAGS",
  "T3_attitude ~ T1_Group * T1_attitude",
  "T3_dehumanization ~ T1_Group * T1_dehumanization",
  "T3_IAT ~ T1_Group * T1_IAT",
  "T3_blastintensity ~ T1_Group * T1_blastintensity",
  "T3_blastduration ~ T1_Group * T1_blastduration",
  "T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration",
  "T3_WHS ~ T1_Group * T1_WHS",
  "T3_CLS ~ T1_Group * T1_CLS"
)

# Make list of all models
models.list2 <- sapply(formulas2, lm, data = data, simplify = FALSE, USE.NAMES = TRUE)

## Attempt with nice_lm_contrasts 
x2 <- nice_lm_contrasts(models.list2, group = "T1_Group", data = data)

x.table2 <- nice_table(x2, highlight = TRUE)
x.table2

Dependent Variable

Comparison

df

t

p

d

95% CI

T3_NOBAGS

Meditation - Waitlist

197

-1.05

.297

-0.24

[-0.57, 0.06]

Reflection - Waitlist

197

0.34

.735

0.05

[-0.31, 0.41]

Meditation - Reflection

197

-1.23

.220

-0.29

[-0.70, 0.09]

T3_attitude

Meditation - Waitlist

197

0.88

.379

0.19

[-0.16, 0.51]

Reflection - Waitlist

197

1.44

.150

0.07

[-0.26, 0.41]

Meditation - Reflection

197

-0.53

.594

0.11

[-0.28, 0.48]

T3_dehumanization

Meditation - Waitlist

196

1.23

.221

0.16

[-0.16, 0.42]

Reflection - Waitlist

196

-0.82

.415

-0.17

[-0.55, 0.22]

Meditation - Reflection

196

1.83

.069

0.33

[-0.05, 0.72]

T3_IAT

Meditation - Waitlist

197

1.30

.196

0.20

[-0.09, 0.55]

Reflection - Waitlist

197

1.15

.250

0.36

[0.02, 0.72]

Meditation - Reflection

197

0.08

.938

-0.16

[-0.55, 0.22]

T3_blastintensity

Meditation - Waitlist

197

-1.79

.076

-0.11

[-0.44, 0.22]

Reflection - Waitlist

197

-2.82

.005**

-0.11

[-0.46, 0.25]

Meditation - Reflection

197

1.02

.310

-0.01

[-0.38, 0.40]

T3_blastduration

Meditation - Waitlist

197

-0.94

.348

-0.06

[-0.36, 0.27]

Reflection - Waitlist

197

-2.69

.008**

-0.21

[-0.58, 0.13]

Meditation - Reflection

197

1.62

.106

0.15

[-0.25, 0.51]

T3_blastintensity.duration

Meditation - Waitlist

197

-1.16

.247

-0.08

[-0.40, 0.23]

Reflection - Waitlist

197

-2.13

.034*

-0.06

[-0.41, 0.32]

Meditation - Reflection

197

0.93

.354

-0.02

[-0.42, 0.35]

T3_WHS

Meditation - Waitlist

197

1.68

.095

0.17

[-0.15, 0.49]

Reflection - Waitlist

197

2.74

.007**

0.11

[-0.23, 0.45]

Meditation - Reflection

197

-1.03

.307

0.06

[-0.30, 0.44]

T3_CLS

Meditation - Waitlist

196

1.99

.048*

0.33

[-0.01, 0.66]

Reflection - Waitlist

196

3.03

.003**

0.37

[-0.02, 0.69]

Meditation - Reflection

196

-0.99

.322

-0.04

[-0.43, 0.33]

Violin Plots & Marginal Means

In this section, we generate various plots: violin plots, marginal means plots, and lighthouse plots, to better illustrate the pairwise contrasts (only for significant contrasts from the previous step).

Time 2

…

Attitudes

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_attitude",
            obs = TRUE,
            signif_annotation = c("*"),
            signif_yposition = 1.5,
            signif_xmin = 2,
            signif_xmax = 3)

The first plot below represents the adjusted (marginal) means, after taking into account covariates such as Time 1 measurements.

means <- estimate_means(models.list$`T2_attitude ~ T1_Group * T1_attitude`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Attitude", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(models.list$`T2_attitude ~ T1_Group * T1_attitude`)
plot(contrasts, estimate_means(models.list$`T2_attitude ~ T1_Group * T1_attitude`)) +
  ylab(paste("Adjusted", "Attitude", "Mean")) +
  theme_modern()

The plot above is a lighthouse plot. These plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

PANAS (Positive Affect)

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_PANAS_pos",
            obs = TRUE,
            comp1 = 1,
            comp2 = 3,
            has.d = TRUE,
            d.x = 1.75,
            d.y = 2)

means <- estimate_means(models.list$`T2_PANAS_pos ~ T1_Group`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Positive Affect", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(models.list$`T2_PANAS_pos ~ T1_Group`)
plot(contrasts, estimate_means(models.list$`T2_PANAS_pos ~ T1_Group`)) +
  ylab(paste("Adjusted", "Positive Affect", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Aggression (Blast Intensity)

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_blastintensity",
            obs = TRUE,
            signif_annotation = c("**", "***"),
            signif_yposition = c(3.5, 4),
            signif_xmin = c(1, 2),
            signif_xmax = c(2, 3))

means <- estimate_means(models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)
plot(contrasts, estimate_means(
  models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)) +
  ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Aggression (Blast Duration)

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_blastduration",
            obs = TRUE,
            signif_annotation = c("**", "***"),
            signif_yposition = c(3.5, 4),
            signif_xmin = c(1, 2),
            signif_xmax = c(2, 3))

means <- estimate_means(models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Duration", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)
plot(contrasts, estimate_means(
  models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)) +
  ylab(paste("Adjusted", "Blast Duration", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Aggression (Blast Intensity * Duration)

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_blastintensity.duration",
            obs = TRUE,
            signif_annotation = c("*", "**"),
            signif_yposition = c(3.5, 4),
            signif_xmin = c(1, 2),
            signif_xmax = c(2, 3))

means <- estimate_means(
  models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
plot(contrasts, estimate_means(
  models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Altruistic Memory

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_memory.altruistic",
            obs = TRUE,
            signif_annotation = c("*"),
            signif_yposition = 3.5,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(
  models.list$`T2_memory.altruistic ~ T1_Group`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "ms to remember altruistic event", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list$`T2_memory.altruistic ~ T1_Group`)
plot(contrasts, estimate_means(
  models.list$`T2_memory.altruistic ~ T1_Group`)) +
  ylab(paste("Adjusted", "ms to remember altruistic event", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Compassion

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_CLS",
            obs = TRUE,
            signif_annotation = c("*", "***"),
            signif_yposition = c(3.2, 2.5),
            signif_xmin = c(1, 2),
            signif_xmax = c(3, 3))

means <- estimate_means(models.list$`T2_CLS ~ T1_Group * T1_CLS`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Compassion", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(models.list$`T2_CLS ~ T1_Group * T1_CLS`)
plot(contrasts, estimate_means(models.list$`T2_CLS ~ T1_Group * T1_CLS`)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Time 3

…

Aggression (Blast Intensity)

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_blastintensity",
            obs = TRUE,
            signif_annotation = c("**"),
            signif_yposition = 2.5,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)
plot(contrasts, estimate_means(
  models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)) +
  ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Aggression (Blast Duration)

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_blastduration",
            obs = TRUE,
            signif_annotation = c("**"),
            signif_yposition = 3,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Duration", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)
plot(contrasts, estimate_means(
  models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)) +
  ylab(paste("Adjusted", "Blast Duration", "Mean")) +
  theme_modern()

Aggression (Blast Intensity * Duration)

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_blastintensity.duration",
            obs = TRUE,
            signif_annotation = c("*"),
            signif_yposition = 3,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(
  models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
plot(contrasts, estimate_means(
  models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Helping

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_WHS",
            obs = TRUE,
            signif_annotation = c("**"),
            signif_yposition = 2.5,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(models.list2$`T3_WHS ~ T1_Group * T1_WHS`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Willingness to Help", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list2$`T3_WHS ~ T1_Group * T1_WHS`)
plot(contrasts, estimate_means(
  models.list2$`T3_WHS ~ T1_Group * T1_WHS`)) +
  ylab(paste("Adjusted", "Willingness to Help", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Compassion

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_CLS",
            obs = TRUE,
            signif_annotation = c("*", "**"),
            signif_yposition = c(3.2, 2.7),
            signif_xmin = c(1, 2),
            signif_xmax = c(3, 3))

means <- estimate_means(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Compassion", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)
plot(contrasts, estimate_means(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).

Compliance Analysis (CACE)

When conducting a Randomized-Controlled Trial (RCT), some participants are randomly assigned to a treatment condition, and others to a control group. However, not everyone assigned to the treatment might follow the treatment protocol (called “treatment compliance”).

According to Sagarin et al. (2014), one sensible approach to address this problem is using the complier average causal effect (CACE), also sometimes known as Local average treatment effect (LATE). According to Wikipedia, it is “the treatment effect for the subset of the sample that takes the treatment if and only if they were assigned to the treatment, otherwise known as the compliers.” In other words, it will be useful if a proportion of your participants assigned to the treatment group did not follow the treatment protocol.

In the following figures, the P > in the compliance column represents the compliance Percentage. The numbers on the right represent Hedge’s g (analogous to Cohen’s d) and its 95% confidence interval.

Time 2

…

Attitudes

data3 <- data2 %>% 
  mutate(part.percent = ifelse(T1_Group == "Waitlist", 0, part.percent),
         part.percent = part.percent * 100,
         T2_memory.altruistic = T2_memory.altruistic * -1,
         T1_GroupReflection = as.numeric(T1_Group == "Reflection"),
         T1_GroupMeditation = as.numeric(T1_Group == "Meditation"),
         across(contains("blast"), \(x) x * -1)) %>% 
  as.data.frame()

caceOutput <- caceSRTBoot(
  T2_attitude ~ T1_attitude + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Attitude (Reflection VS Waitlist)")

PANAS (Positive Affect)

caceOutput <- caceSRTBoot(
  T2_PANAS_pos ~ T1_GroupMeditation,
  intervention = "T1_GroupMeditation",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Positive Affect (Meditation VS Waitlist)")

Aggression (Blast Intensity)

caceOutput <- caceSRTBoot(
  T2_blastintensity ~ T1_blastintensity + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Intensity (Reflection VS Waitlist)")

Aggression (Blast Duration)

caceOutput <- caceSRTBoot(
  T2_blastduration ~ T1_blastduration + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Duration (Reflection VS Waitlist)")

Aggression (Blast Intensity * Duration)

caceOutput <- caceSRTBoot(
  T2_blastintensity.duration ~ T1_blastintensity.duration + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Intensity * Duration (Reflection VS Waitlist)")

Altruistic Memory

caceOutput <- caceSRTBoot(
  T2_memory.altruistic ~ T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "ms to remember altruistic event (Reflection VS Waitlist)")

Compassion

caceOutput <- caceSRTBoot(
  T2_CLS ~ T1_CLS + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T2_CLS ~ T1_CLS + T1_GroupMeditation,
  intervention = "T1_GroupMeditation",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Meditation VS Waitlist)")

Time 3

…

Attitudes

caceOutput <- caceSRTBoot(
  T3_attitude ~ T1_attitude + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Attitude (Reflection VS Waitlist)")

Aggression (Blast Intensity)

caceOutput <- caceSRTBoot(
  T3_blastintensity ~ T1_blastintensity + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Intensity (Reflection VS Waitlist)")

Aggression (Blast Duration)

caceOutput <- caceSRTBoot(
  T3_blastduration ~ T1_blastduration + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Duration (Reflection VS Waitlist)")

Aggression (Blast Intensity * Duration)

caceOutput <- caceSRTBoot(
  T3_blastintensity.duration ~ T1_blastintensity.duration + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Intensity * Duration (Reflection VS Waitlist)")

Helping

caceOutput <- caceSRTBoot(
  T3_WHS ~ T1_WHS + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Willingness to Help (Reflection VS Waitlist)")

Compassion

caceOutput <- caceSRTBoot(
  T3_CLS ~ T1_CLS + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T3_CLS ~ T1_CLS + T1_GroupMeditation,
  intervention = "T1_GroupMeditation",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Meditation VS Waitlist)")

What we see is general trend toward larger effects with higher compliance, although this trend does not appear exceptionally strong.

Moderations (Ego Depletion)

…

It seems like only IAT (implicit aggression) and not NOBAGS (explicit attitude toward aggression) is a significant moderator.

Moderator: IAT*

blastintensity.duration

# CREATE OUR DUMMY VARIABLES FOR T1_Group!

data$T1_GroupReflection <- as.numeric(data$T1_Group == "Reflection")
data$T1_GroupMeditation <- as.numeric(data$T1_Group == "Meditation")

################################################
################################################

# T2_blastintensity.duration - IAT
T2_blastintensity.duration.IAT <- lm(
  T2_blastintensity.duration ~ T1_GroupReflection + T1_GroupMeditation +
    T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + T2_IAT:T2_Condition + 
    T1_GroupReflection:T2_IAT:T2_Condition + T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_blastintensity.duration.IAT)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_blastintensity.duration.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_blastintensity.duration

T1_GroupReflection

191

-405.54

-0.31

.760

.00

[0.00, 0.01]

T1_GroupMeditation

191

-54.71

-0.03

.974

.00

[0.00, 0.00]

T2_IAT

191

-155.01

-0.11

.915

.00

[0.00, 0.00]

T2_ConditionDepleted

191

2,409.51

1.68

.096

.01

[0.00, 0.04]

T1_GroupReflection × T2_IAT

191

215.91

0.10

.917

.00

[0.00, 0.00]

T1_GroupMeditation × T2_IAT

191

-69.14

-0.03

.979

.00

[0.00, 0.00]

T1_GroupReflection × T2_ConditionDepleted

191

-2,181.85

-0.85

.398

.00

[0.00, 0.02]

T1_GroupMeditation × T2_ConditionDepleted

191

-5,380.99

-2.14

.033*

.02

[0.00, 0.06]

T2_IAT × T2_ConditionDepleted

191

2,746.44

1.41

.161

.01

[0.00, 0.04]

T1_GroupReflection × T2_IAT × T2_ConditionDepleted

191

-3,823.21

-1.06

.290

.01

[0.00, 0.03]

T1_GroupMeditation × T2_IAT × T2_ConditionDepleted

191

-7,357.76

-2.06

.041*

.02

[0.00, 0.06]

# Make interaction plot
interact_plot(T2_blastintensity.duration.IAT, pred = "T2_IAT", modx = "T2_Condition",
              interval = TRUE, mod2 = "T1_GroupMeditation",
              x.label = "Implicit Aggression (IAT)",
              y.label = "Blast Intensity * Duration (Taylor)", 
              mod2.labels = c("Waitlist Group", "Meditation Group"),
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

What we see here is that for the waitlist group, there is no interaction between being depleted and implicit aggression, whereas for the meditation group, people become less aggressive after being depleted, the higher their implicit aggression.

blastintensity

# T2_blastintensity - IAT
T2_blastintensity.IAT <- lm(
  T2_blastintensity ~ T1_GroupReflection + T1_GroupMeditation +
    T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + 
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_blastintensity.IAT)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_blastintensity.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_blastintensity

T1_GroupReflection

191

-0.59

-0.67

.502

.00

[0.00, 0.02]

T1_GroupMeditation

191

-0.19

-0.17

.864

.00

[0.00, 0.00]

T2_IAT

191

-0.09

-0.09

.928

.00

[0.00, 0.00]

T2_ConditionDepleted

191

0.84

0.88

.381

.00

[0.00, 0.02]

T1_GroupReflection × T2_IAT

191

-0.07

-0.05

.961

.00

[0.00, 0.00]

T1_GroupMeditation × T2_IAT

191

-0.07

-0.04

.966

.00

[0.00, 0.00]

T1_GroupReflection × T2_ConditionDepleted

191

-0.79

-0.46

.646

.00

[0.00, 0.01]

T1_GroupMeditation × T2_ConditionDepleted

191

-2.29

-1.38

.171

.01

[0.00, 0.04]

T2_IAT × T2_ConditionDepleted

191

1.08

0.83

.406

.00

[0.00, 0.02]

T1_GroupReflection × T2_IAT × T2_ConditionDepleted

191

-1.88

-0.79

.432

.00

[0.00, 0.02]

T1_GroupMeditation × T2_IAT × T2_ConditionDepleted

191

-3.49

-1.47

.142

.01

[0.00, 0.04]

blastduration *

# T2_blastduration - IAT
T2_blastduration.IAT <- lm(
  T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
    T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_blastduration.IAT)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_blastduration.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_blastduration

T1_GroupReflection

191

-151.94

-0.84

.403

.00

[0.00, 0.02]

T1_GroupMeditation

191

133.46

0.58

.560

.00

[0.00, 0.01]

T2_IAT

191

-92.44

-0.47

.641

.00

[0.00, 0.01]

T2_ConditionDepleted

191

271.17

1.38

.169

.01

[0.00, 0.04]

T1_GroupReflection × T2_IAT

191

8.01

0.03

.977

.00

[0.00, 0.00]

T1_GroupMeditation × T2_IAT

191

196.17

0.55

.582

.00

[0.00, 0.01]

T1_GroupReflection × T2_ConditionDepleted

191

-231.85

-0.66

.511

.00

[0.00, 0.01]

T1_GroupMeditation × T2_ConditionDepleted

191

-706.58

-2.06

.041*

.02

[0.00, 0.06]

T2_IAT × T2_ConditionDepleted

191

326.41

1.22

.223

.01

[0.00, 0.03]

T1_GroupReflection × T2_IAT × T2_ConditionDepleted

191

-471.58

-0.96

.339

.00

[0.00, 0.02]

T1_GroupMeditation × T2_IAT × T2_ConditionDepleted

191

-962.19

-1.97

.050

.02

[0.00, 0.06]

# Make interaction plot
interact_plot(T2_blastduration.IAT, pred = "T2_IAT", modx = "T2_Condition",
              interval = TRUE, mod2 = "T1_GroupMeditation",
              x.label = "Implicit Aggression (IAT)",
              y.label = "Blast Duration (Taylor)", 
              mod2.labels = c("Waitlist Group", "Meditation Group"),
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

What we see here is that for the waitlist group, there is no interaction between being depleted and implicit aggression, whereas for the meditation group, people become less aggressive after being depleted, the higher their implicit aggression. However, that variable (blastduration alone) was not in the preregistration, so we might not report this finding.

Charity

# T2_Charity - IAT
T2_Charity.IAT <- lm(
  T2_Charity ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
    T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +  
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition + T2_Familiarity, 
  data = data)

check_model(T2_Charity.IAT)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_Charity.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_Charity

T1_GroupReflection

182

0.91

0.71

.480

.00

[0.00, 0.02]

T1_GroupMeditation

182

-0.72

-0.43

.670

.00

[0.00, 0.01]

T2_IAT

182

-0.47

-0.34

.732

.00

[0.00, 0.01]

T2_ConditionDepleted

182

1.04

0.76

.451

.00

[0.00, 0.02]

T2_Familiarity

182

1.04

3.56

< .001***

.06

[0.00, 0.13]

T1_GroupReflection × T2_IAT

182

1.55

0.78

.435

.00

[0.00, 0.02]

T1_GroupMeditation × T2_IAT

182

-0.26

-0.10

.921

.00

[0.00, 0.00]

T1_GroupReflection × T2_ConditionDepleted

182

-2.03

-0.79

.430

.00

[0.00, 0.02]

T1_GroupMeditation × T2_ConditionDepleted

182

-0.31

-0.13

.900

.00

[0.00, 0.00]

T2_IAT × T2_ConditionDepleted

182

1.32

0.70

.487

.00

[0.00, 0.02]

T1_GroupReflection × T2_IAT × T2_ConditionDepleted

182

-0.78

-0.22

.825

.00

[0.00, 0.00]

T1_GroupMeditation × T2_IAT × T2_ConditionDepleted

182

0.12

0.03

.974

.00

[0.00, 0.00]

Compassion*

# Compassion - IAT
T2_CLS.IAT <- lm(
  T2_CLS ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + 
    T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_CLS.IAT)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_CLS.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_CLS

T1_GroupReflection

191

1.43

3.77

< .001***

.06

[0.00, 0.12]

T1_GroupMeditation

191

1.42

2.97

.003**

.04

[0.00, 0.09]

T2_IAT

191

-1.47

-3.54

< .001***

.06

[0.00, 0.11]

T2_ConditionDepleted

191

1.15

2.80

.006**

.03

[0.00, 0.08]

T1_GroupReflection × T2_IAT

191

1.69

2.86

.005**

.04

[0.00, 0.08]

T1_GroupMeditation × T2_IAT

191

1.77

2.37

.019*

.02

[0.00, 0.06]

T1_GroupReflection × T2_ConditionDepleted

191

0.09

0.12

.907

.00

[0.00, 0.00]

T1_GroupMeditation × T2_ConditionDepleted

191

-1.87

-2.61

.010**

.03

[0.00, 0.07]

T2_IAT × T2_ConditionDepleted

191

1.96

3.51

.001***

.05

[0.00, 0.11]

T1_GroupReflection × T2_IAT × T2_ConditionDepleted

191

-0.03

-0.03

.979

.00

[0.00, 0.00]

T1_GroupMeditation × T2_IAT × T2_ConditionDepleted

191

-3.07

-3.00

.003**

.04

[0.00, 0.09]

# Make interaction plot
interact_plot(T2_CLS.IAT, pred = "T2_IAT", modx = "T2_Condition",
              interval = TRUE, mod2 = "T1_GroupMeditation",
              x.label = "Implicit Aggression (IAT)",
              y.label = "Compassionate Love",
              mod2.labels = c("Waitlist Group", "Meditation Group"),
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

What we see here is that for the waitlist group, the effect of implicit aggression clearly depends on depletion: implicit aggression leads to lower compassion in the control condition (expected), but to higher compassion in the depletion group (unexpected). Whereas for the meditation group, it seems the depletion has little effect on compassion (but the trend is reversed relative to the waitlist group).

Helping

# T2_WHS - IAT
T2_WHS.IAT <- lm(
  T2_WHS ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + T2_Condition + 
    T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT + T1_GroupReflection:T2_Condition +
    T1_GroupMeditation:T2_Condition + T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_WHS.IAT)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_WHS.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_WHS

T1_GroupReflection

191

0.53

1.43

.154

.01

[0.00, 0.04]

T1_GroupMeditation

191

0.69

1.48

.141

.01

[0.00, 0.04]

T2_IAT

191

-0.73

-1.81

.072

.02

[0.00, 0.05]

T2_ConditionDepleted

191

0.50

1.25

.213

.01

[0.00, 0.03]

T1_GroupReflection × T2_IAT

191

0.90

1.57

.119

.01

[0.00, 0.04]

T1_GroupMeditation × T2_IAT

191

0.95

1.31

.193

.01

[0.00, 0.03]

T1_GroupReflection × T2_ConditionDepleted

191

0.38

0.53

.595

.00

[0.00, 0.01]

T1_GroupMeditation × T2_ConditionDepleted

191

-0.88

-1.26

.209

.01

[0.00, 0.03]

T2_IAT × T2_ConditionDepleted

191

0.85

1.55

.123

.01

[0.00, 0.04]

T1_GroupReflection × T2_IAT × T2_ConditionDepleted

191

0.49

0.49

.624

.00

[0.00, 0.01]

T1_GroupMeditation × T2_IAT × T2_ConditionDepleted

191

-1.25

-1.25

.212

.01

[0.00, 0.03]

Altruistic Memory*

# T2_memory.altruistic - IAT
T2_memory.altruistic.IAT <- lm(
  T2_memory.altruistic ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + 
    T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +  
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_memory.altruistic.IAT)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_memory.altruistic.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_memory.altruistic

T1_GroupReflection

191

-3,870.57

-2.08

.039*

.02

[0.00, 0.06]

T1_GroupMeditation

191

-4,129.04

-1.76

.081

.02

[0.00, 0.05]

T2_IAT

191

2,529.17

1.24

.216

.01

[0.00, 0.03]

T2_ConditionDepleted

191

-4,289.34

-2.12

.035*

.02

[0.00, 0.06]

T1_GroupReflection × T2_IAT

191

-2,907.13

-1.00

.318

.00

[0.00, 0.02]

T1_GroupMeditation × T2_IAT

191

-5,638.28

-1.54

.126

.01

[0.00, 0.04]

T1_GroupReflection × T2_ConditionDepleted

191

3,963.65

1.10

.275

.01

[0.00, 0.03]

T1_GroupMeditation × T2_ConditionDepleted

191

4,949.85

1.40

.162

.01

[0.00, 0.04]

T2_IAT × T2_ConditionDepleted

191

-5,349.23

-1.95

.053

.02

[0.00, 0.05]

T1_GroupReflection × T2_IAT × T2_ConditionDepleted

191

5,764.69

1.14

.256

.01

[0.00, 0.03]

T1_GroupMeditation × T2_IAT × T2_ConditionDepleted

191

6,772.42

1.35

.179

.01

[0.00, 0.03]

# Make interaction plot
interact_plot(T2_memory.altruistic.IAT, pred = "T2_IAT", modx = "T2_Condition",
              interval = TRUE, mod2 = "T1_GroupMeditation",
              x.label = "Implicit Aggression (IAT)",
              y.label = "Reaction Time (Altruistic Memory)",
              mod2.labels = c("Waitlist Group", "Meditation Group"),
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

What we see here is that for the waitlist group, the depletion seems to have little effect on reaction time to remember an altruistic event. Whereas for the meditation group, higher implicit aggression relates to shorter reaction time (unexpected), unless they are depleted (but the trend is reversed relative to the waitlist group). This result is somewhat unexpected, but it could be that in their case, the heart takes over the mind.

Aggressive Memory

# T2_memory.aggressive - IAT
T2_memory.aggressive.IAT <- lm(
  T2_memory.aggressive ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + 
    T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +  
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_memory.aggressive.IAT)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_memory.aggressive.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_memory.aggressive

T1_GroupReflection

191

2,692.41

0.61

.545

.00

[0.00, 0.01]

T1_GroupMeditation

191

4,552.33

0.81

.417

.00

[0.00, 0.02]

T2_IAT

191

-2,085.63

-0.43

.668

.00

[0.00, 0.01]

T2_ConditionDepleted

191

1,809.33

0.38

.707

.00

[0.00, 0.01]

T1_GroupReflection × T2_IAT

191

1,263.35

0.18

.855

.00

[0.00, 0.00]

T1_GroupMeditation × T2_IAT

191

5,654.74

0.65

.518

.00

[0.00, 0.01]

T1_GroupReflection × T2_ConditionDepleted

191

-2,341.73

-0.27

.786

.00

[0.00, 0.01]

T1_GroupMeditation × T2_ConditionDepleted

191

-7,225.39

-0.86

.391

.00

[0.00, 0.02]

T2_IAT × T2_ConditionDepleted

191

-87.36

-0.01

.989

.00

[0.00, 0.00]

T1_GroupReflection × T2_IAT × T2_ConditionDepleted

191

224.90

0.02

.985

.00

[0.00, 0.00]

T1_GroupMeditation × T2_IAT × T2_ConditionDepleted

191

-14,323.47

-1.20

.233

.01

[0.00, 0.03]

Moderator: NOBAGS

blastintensity.duration

# T2_blastintensity.duration - NOBAGS
T2_blastintensity.duration.NOBAGS <- lm(
  T2_blastintensity.duration ~ T1_GroupReflection + T1_GroupMeditation +
    T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + 
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_blastintensity.duration.NOBAGS)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_blastintensity.duration.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_blastintensity.duration

T1_GroupReflection

191

-3,581.79

-1.16

.246

.01

[0.00, 0.03]

T1_GroupMeditation

191

-974.72

-0.29

.770

.00

[0.00, 0.01]

T2_NOBAGS

191

-733.10

-0.70

.488

.00

[0.00, 0.02]

T2_ConditionDepleted

191

-5,165.29

-1.60

.112

.01

[0.00, 0.04]

T1_GroupReflection × T2_NOBAGS

191

1,525.17

1.02

.309

.01

[0.00, 0.02]

T1_GroupMeditation × T2_NOBAGS

191

451.50

0.28

.781

.00

[0.00, 0.01]

T1_GroupReflection × T2_ConditionDepleted

191

8,010.38

1.63

.104

.01

[0.00, 0.04]

T1_GroupMeditation × T2_ConditionDepleted

191

1,410.90

0.29

.769

.00

[0.00, 0.01]

T2_NOBAGS × T2_ConditionDepleted

191

2,862.52

1.82

.070

.02

[0.00, 0.05]

T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted

191

-3,841.98

-1.56

.121

.01

[0.00, 0.04]

T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted

191

-613.81

-0.25

.800

.00

[0.00, 0.01]

blastintensity

# T2_blastintensity - NOBAGS
T2_blastintensity.NOBAGS <- lm(
  T2_blastintensity ~ T1_GroupReflection + T1_GroupMeditation +
    T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + 
    T1_GroupMeditation:T2_NOBAGS +  T1_GroupReflection:T2_Condition + 
    T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition + 
    T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_blastintensity.NOBAGS)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_blastintensity.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_blastintensity

T1_GroupReflection

191

-2.72

-1.35

.177

.01

[0.00, 0.03]

T1_GroupMeditation

191

-0.47

-0.21

.831

.00

[0.00, 0.00]

T2_NOBAGS

191

-0.55

-0.80

.422

.00

[0.00, 0.02]

T2_ConditionDepleted

191

-3.82

-1.81

.072

.02

[0.00, 0.05]

T1_GroupReflection × T2_NOBAGS

191

1.07

1.10

.274

.01

[0.00, 0.03]

T1_GroupMeditation × T2_NOBAGS

191

0.14

0.13

.899

.00

[0.00, 0.00]

T1_GroupReflection × T2_ConditionDepleted

191

4.54

1.42

.158

.01

[0.00, 0.04]

T1_GroupMeditation × T2_ConditionDepleted

191

0.70

0.22

.824

.00

[0.00, 0.00]

T2_NOBAGS × T2_ConditionDepleted

191

1.96

1.91

.057

.02

[0.00, 0.05]

T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted

191

-1.97

-1.22

.223

.01

[0.00, 0.03]

T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted

191

-0.11

-0.07

.947

.00

[0.00, 0.00]

blastduration*

# T2_blastduration - NOBAGS
T2_blastduration.NOBAGS <- lm(
  T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
    T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + 
    T1_GroupMeditation:T2_NOBAGS +  T1_GroupReflection:T2_Condition +
    T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition +
    T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_blastduration.NOBAGS)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_blastduration.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_blastduration

T1_GroupReflection

191

-491.94

-1.20

.233

.01

[0.00, 0.03]

T1_GroupMeditation

191

26.08

0.06

.953

.00

[0.00, 0.00]

T2_NOBAGS

191

-119.38

-0.85

.397

.00

[0.00, 0.02]

T2_ConditionDepleted

191

-822.26

-1.90

.058

.02

[0.00, 0.05]

T1_GroupReflection × T2_NOBAGS

191

157.43

0.79

.431

.00

[0.00, 0.02]

T1_GroupMeditation × T2_NOBAGS

191

-8.25

-0.04

.970

.00

[0.00, 0.00]

T1_GroupReflection × T2_ConditionDepleted

191

738.07

1.13

.261

.01

[0.00, 0.03]

T1_GroupMeditation × T2_ConditionDepleted

191

-243.02

-0.38

.705

.00

[0.00, 0.01]

T2_NOBAGS × T2_ConditionDepleted

191

438.08

2.09

.038*

.02

[0.00, 0.06]

T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted

191

-302.62

-0.92

.359

.00

[0.00, 0.02]

T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted

191

143.03

0.44

.659

.00

[0.00, 0.01]

# Make interaction plot
interact_plot(T2_blastduration.NOBAGS, pred = "T2_NOBAGS", modx = "T2_Condition",
              interval = TRUE, 
              x.label = "Normative beliefs about aggression (NOBAGS)",
              y.label = "Blast Duration (Taylor)",
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

What we see here is that there seems to be no relation between attitude toward aggression on behavioural aggression in the control condition, whereas when depleted, attitude toward aggression relates to higher behavioural aggression. Although this result is theoretically consistent with the literature, it is likely a false positive given our high number of tests, the fact that this is the only variable that NOBAGS moderates, and that the p value is relatively close to 0.5. Furthermore that variable (blastduration alone) was not in the preregistration, so we might not report this finding.

Charity

# T2_Charity - NOBAGS
T2_Charity.NOBAGS <- lm(
  T2_Charity ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS +
    T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition + T2_Familiarity, 
  data = data)

check_model(T2_Charity.NOBAGS)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_Charity.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_Charity

T1_GroupReflection

182

1.16

0.39

.695

.00

[0.00, 0.01]

T1_GroupMeditation

182

1.13

0.36

.721

.00

[0.00, 0.01]

T2_NOBAGS

182

0.62

0.62

.534

.00

[0.00, 0.01]

T2_ConditionDepleted

182

2.98

0.96

.339

.00

[0.00, 0.02]

T2_Familiarity

182

1.17

3.93

< .001***

.08

[0.01, 0.15]

T1_GroupReflection × T2_NOBAGS

182

-0.51

-0.35

.725

.00

[0.00, 0.01]

T1_GroupMeditation × T2_NOBAGS

182

-0.86

-0.56

.577

.00

[0.00, 0.01]

T1_GroupReflection × T2_ConditionDepleted

182

-8.02

-1.69

.092

.01

[0.00, 0.05]

T1_GroupMeditation × T2_ConditionDepleted

182

-5.03

-1.09

.275

.01

[0.00, 0.03]

T2_NOBAGS × T2_ConditionDepleted

182

-1.38

-0.92

.359

.00

[0.00, 0.02]

T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted

182

3.29

1.37

.174

.01

[0.00, 0.03]

T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted

182

2.40

1.03

.303

.01

[0.00, 0.02]

Compassion

# Compassion - NOBAGS
T2_CLS.NOBAGS <- lm(
  T2_CLS ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS + 
    T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +  
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_CLS.NOBAGS)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_CLS.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_CLS

T1_GroupReflection

191

0.13

0.14

.889

.00

[0.00, 0.00]

T1_GroupMeditation

191

1.51

1.54

.125

.01

[0.00, 0.04]

T2_NOBAGS

191

-0.24

-0.76

.448

.00

[0.00, 0.02]

T2_ConditionDepleted

191

0.63

0.67

.505

.00

[0.00, 0.01]

T1_GroupReflection × T2_NOBAGS

191

0.14

0.31

.757

.00

[0.00, 0.01]

T1_GroupMeditation × T2_NOBAGS

191

-0.61

-1.27

.204

.01

[0.00, 0.03]

T1_GroupReflection × T2_ConditionDepleted

191

-0.27

-0.19

.850

.00

[0.00, 0.00]

T1_GroupMeditation × T2_ConditionDepleted

191

-1.68

-1.19

.235

.01

[0.00, 0.03]

T2_NOBAGS × T2_ConditionDepleted

191

-0.38

-0.82

.414

.00

[0.00, 0.02]

T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted

191

0.10

0.14

.892

.00

[0.00, 0.00]

T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted

191

0.88

1.24

.215

.01

[0.00, 0.03]

Helping

# T2_WHS - NOBAGS
T2_WHS.NOBAGS <- lm(
  T2_WHS ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS + T2_Condition + 
    T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS + 
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + 
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_WHS.NOBAGS)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_WHS.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_WHS

T1_GroupReflection

191

-1.04

-1.22

.222

.01

[0.00, 0.03]

T1_GroupMeditation

191

0.92

1.01

.316

.00

[0.00, 0.02]

T2_NOBAGS

191

-0.47

-1.63

.104

.01

[0.00, 0.04]

T2_ConditionDepleted

191

-1.67

-1.87

.063

.02

[0.00, 0.05]

T1_GroupReflection × T2_NOBAGS

191

0.51

1.23

.219

.01

[0.00, 0.03]

T1_GroupMeditation × T2_NOBAGS

191

-0.43

-0.97

.335

.00

[0.00, 0.02]

T1_GroupReflection × T2_ConditionDepleted

191

1.90

1.41

.161

.01

[0.00, 0.04]

T1_GroupMeditation × T2_ConditionDepleted

191

-0.54

-0.41

.683

.00

[0.00, 0.01]

T2_NOBAGS × T2_ConditionDepleted

191

0.81

1.87

.063

.02

[0.00, 0.05]

T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted

191

-0.97

-1.43

.154

.01

[0.00, 0.04]

T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted

191

0.26

0.40

.693

.00

[0.00, 0.01]

Altruistic Memory

# T2_memory.altruistic - NOBAGS
T2_memory.altruistic.NOBAGS <- lm(
  T2_memory.altruistic ~ T1_GroupReflection + T1_GroupMeditation +
    T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + 
    T1_GroupMeditation:T2_NOBAGS + T1_GroupReflection:T2_Condition + 
    T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition + 
    T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_memory.altruistic.NOBAGS)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_memory.altruistic.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_memory.altruistic

T1_GroupReflection

191

-1,307.87

-0.30

.761

.00

[0.00, 0.01]

T1_GroupMeditation

191

1,767.77

0.38

.705

.00

[0.00, 0.01]

T2_NOBAGS

191

1,260.19

0.86

.393

.00

[0.00, 0.02]

T2_ConditionDepleted

191

-3,774.82

-0.84

.404

.00

[0.00, 0.02]

T1_GroupReflection × T2_NOBAGS

191

-360.53

-0.17

.863

.00

[0.00, 0.00]

T1_GroupMeditation × T2_NOBAGS

191

-1,283.78

-0.57

.572

.00

[0.00, 0.01]

T1_GroupReflection × T2_ConditionDepleted

191

7,532.92

1.10

.273

.01

[0.00, 0.03]

T1_GroupMeditation × T2_ConditionDepleted

191

1,944.98

0.29

.772

.00

[0.00, 0.01]

T2_NOBAGS × T2_ConditionDepleted

191

1,527.83

0.70

.487

.00

[0.00, 0.02]

T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted

191

-3,795.29

-1.10

.272

.01

[0.00, 0.03]

T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted

191

-415.11

-0.12

.903

.00

[0.00, 0.00]

Aggressive Memory

# T2_memory.aggressive - NOBAGS
T2_memory.aggressive.NOBAGS <- lm(
  T2_memory.aggressive ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS + 
    T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + 
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_memory.aggressive.NOBAGS)
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

T2_memory.aggressive.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

Dependent Variable

Predictor

df

b

t

p

sr2

95% CI

T2_memory.aggressive

T1_GroupReflection

191

8,077.14

0.79

.429

.00

[0.00, 0.02]

T1_GroupMeditation

191

-5,022.93

-0.46

.650

.00

[0.00, 0.01]

T2_NOBAGS

191

-1,512.22

-0.43

.665

.00

[0.00, 0.01]

T2_ConditionDepleted

191

1,092.80

0.10

.919

.00

[0.00, 0.00]

T1_GroupReflection × T2_NOBAGS

191

-3,290.98

-0.67

.506

.00

[0.00, 0.01]

T1_GroupMeditation × T2_NOBAGS

191

3,160.31

0.59

.557

.00

[0.00, 0.01]

T1_GroupReflection × T2_ConditionDepleted

191

-17,673.12

-1.09

.278

.01

[0.00, 0.03]

T1_GroupMeditation × T2_ConditionDepleted

191

-4,542.27

-0.29

.775

.00

[0.00, 0.01]

T2_NOBAGS × T2_ConditionDepleted

191

395.27

0.08

.939

.00

[0.00, 0.00]

T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted

191

8,078.65

0.99

.323

.00

[0.00, 0.02]

T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted

191

4,027.85

0.50

.616

.00

[0.00, 0.01]

Attitudes & Charity

…

Let us explore some differences in attitudes toward various social groups and charities.

Social Groups

data %>%
  select(starts_with("T1_attitude")) %>% 
  get_label %>% 
  lapply(function(x) gsub(".*- ", "", x)) %>% 
  unlist() %>% unname
##  [1] "Personnes noires"    "Personnes sans-abri" "Autochtones"        
##  [4] "Musulmans"           "Réfugiés"            "Femmes"             
##  [7] "Animaux"             "Personnes âgées"     "Personnes blanches" 
## [10] ""
social.groups <- c("Blacks", "Homeless", "Native", "Muslims", "Refugees", "Women",
                   "Animals", "Elderly", "Whites")

charities <- data %>%
  select(ends_with("1_1") & contains("charity")) %>% 
  get_label %>% 
  lapply(function(x) gsub(".*- ", "", x)) %>% 
  unlist() %>% unname
charities
##  [1] "Moisson Montréal"                                                
##  [2] "Accueil Bonneau"                                                 
##  [3] "Association de Montréal pour la déficience intellectuelle (AMDI)"
##  [4] "Centre d'amitié autochtone de Montréal (CAAM)"                   
##  [5] "Fondation du Centre des femmes de Montréal"                      
##  [6] "Centre communautaire LGBTQ+ de Montréal"                         
##  [7] "Équiterre"                                                       
##  [8] "Club des petits déjeuners du Québec"                             
##  [9] "Fondation québécoise du cancer"                                  
## [10] "Société Parkinson du Québec"                                     
## [11] "Oxfam-Québec"                                                    
## [12] "Fondation Papillon"                                              
## [13] "Croix-Rouge canadienne"                                          
## [14] "Centraide"                                                       
## [15] "Société canadienne de l'autisme"                                 
## [16] "Fondation David Suzuki"                                          
## [17] "Conservation de la nature Canada"                                
## [18] "Fondation des maladies du cœur et de l'AVC du Canada"            
## [19] "UNICEF"                                                          
## [20] "Amnistie internationale"                                         
## [21] "Greenpeace"                                                      
## [22] "Fonds mondial pour la nature (WWF)"                              
## [23] "Médecins sans frontières (MSF)"                                  
## [24] "Armée du Salut"
regions <- c("Montreal", "Quebec", "Canada", "International")
regions
## [1] "Montreal"      "Quebec"        "Canada"        "International"

Time 1

Explicit Attitude

data %>%
  select(T1_attitude_1:T1_attitude_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Positive Explicit Attitude",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

Dehumanization

data %>%
  select(T1_dehumanization_1:T1_dehumanization_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Humanization",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

Time 2

Explicit Attitude

data %>%
  select(T2_attitude_1:T2_attitude_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Positive Explicit Attitude",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

Dehumanization

data %>%
  select(T2_dehumanization_1:T2_dehumanization_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Humanization",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

Amount Donated

data %>%
  select(contains("charity") & ends_with("1_1")) %>% 
  pivot_longer(cols = everything(), 
               names_to = "charity",
               values_to = "donation") %>% 
  mutate(charity = factor(charity, labels = charities)) %>% 
  nice_violin(group = "charity",
              response = "donation",
              ytitle = "Amount Donated",
              CIcap.width = 0.5,
              obs = "jitter",
              order.groups = "increasing",
              xlabels.angle = 75)

Charity Familiarity

data %>%
  select(contains("charity") & ends_with("2_1")) %>% 
  pivot_longer(cols = everything(), 
               names_to = "charity",
               values_to = "familiarity") %>% 
  mutate(charity = factor(charity, labels = charities)) %>% 
  nice_violin(group = "charity",
              response = "familiarity",
              ytitle = "Familiarity with Charity",
              CIcap.width = 0.5,
              obs = "jitter",
              order.groups = "increasing",
              xlabels.angle = 75)

Familiarity & Donation

data %>%
  nice_scatter(predictor = "T2_Familiarity",
               response = "T2_Charity",
               ytitle = "Donation Amount",
               xtitle = "Familiarity with Charity",
               has.jitter = TRUE,
               has.legend = TRUE,
               has.r = TRUE,
               has.p = TRUE)

Region & Donation

data %>%
  select(contains("charity") & ends_with("1_1")) %>% 
  pivot_longer(cols = everything(), 
               names_to = "charity",
               values_to = "donation") %>% 
  mutate(charity = factor(charity, labels = charities),
         region = case_match(
           charity,
           charities[1:6] ~ regions[1],
           charities[7:12] ~ regions[2],
           charities[13:18] ~ regions[3],
           charities[19:24] ~ regions[4]
         )) %>% 
  nice_violin(group = "region",
              response = "donation",
              ytitle = "Amount Donated",
              obs = "jitter",
              order.groups = "increasing")

Region & Familiarity

data %>%
  select(contains("charity") & ends_with("2_1")) %>% 
  pivot_longer(cols = everything(), 
               names_to = "charity",
               values_to = "familiarity") %>% 
  mutate(charity = factor(charity, labels = charities),
         region = case_match(
           charity,
           charities[1:6] ~ regions[1],
           charities[7:12] ~ regions[2],
           charities[13:18] ~ regions[3],
           charities[19:24] ~ regions[4]
         )) %>% 
  nice_violin(group = "region",
              response = "familiarity",
              ytitle = "Familiarity with Charity",
              obs = "jitter",
              order.groups = "increasing")

Time 3

Explicit Attitude

data %>%
  select(T3_attitude_1:T3_attitude_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Positive Explicit Attitude",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

Dehumanization

data %>%
  select(T3_dehumanization_1:T3_dehumanization_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Humanization",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

Discussion

In this report, we aimed to compare two types of loving-kindness meditation, one more embodied, based on meditation, and one more cognitive, based on intellectual reflection, to a waitlist control group. We compared those groups on several variables relating to prosociality (i.e., on affect, attitude and behaviour). Groups were measured three times: Week 0 (T1), Week 6 (T2), and Week 13 (T3) so it was able to compare for baseline but also see how robust the effects, if any, are through times. We were also interested in assessing whether these effects depend on other personality variables (i.e., moderators).

Group Differences at Time 2

Our contrasts analyses first revealed group differences at Time 2. Both the meditation and reflection groups showed moderately more compassionate love than the waitlist group, but only the meditation group showed moderately more positive affect than the waitlist group. However, the reflection group showed a little more positive explicit attitudes toward various social groups, as well as moderately shorter reaction times to remember an altruistic event than the waitlist group (suggesting that altruism was more cognitively accessible to them). Furthermore, the reflection group showed a little lower behavioural aggression (blast intensity, blast duration, and blast intensity * duration) than both the waitlist group and the meditation group.

Group Differences at Time 3

Our contrasts analyses also revealed group differences at Time 3. However, only the reflection group showed lasting positive effects on attitudes (still small effect), behavioural aggression (still small effect), and compassion (still moderate effect), suggesting these effects are durable in time. Furthermore, the reflection group showed a delayed onset effect on willingness to help, whereas they were a little more willing to help in various hypothetical scenarios than the control group.

Moderations at Time 2

First, attitudes toward aggression (NOBAGS) only moderated one variable, blast duration. In short, while NOBAGS does not affect blast duration in the control condition, it relates to higher blast duration in the depletion condition. Although this result is theoretically consistent with the literature, it is likely a false positive given our high number of tests, the fact that this is the only variable that NOBAGS moderates, and that the p value is relatively close to 0.5. Furthermore that variable (blastduration alone) was not in the preregistration, so we might not report this finding.

Second, implicit aggression (IAT) seems to have moderated several variables. Like for NOBAGS, it also moderated blast duration, but in a three-way interaction this time. Surprisingly, for the meditation group, implicit aggression related to lower aggression, but only when depleted, whereas there was no such interaction in the waitlist group. However, as mentioned before, that variable (blastduration alone) was not in the preregistration, so we might not report this finding.

Third, implicit aggression also moderated compassionate love, again in a three-way interaction. For the waitlist group, the effect of implicit aggression clearly depends on depletion: implicit aggression relates to lower compassion in the control group (expected), but to higher compassion in the depletion group (unexpected). However, for the meditation group, the effect was absent or partly reversed.

Fourth, implicit aggression also moderated reaction time to remember an altruistic event, again in a three-way interaction. For the meditation group, higher implicit aggression relates to shorter reaction time (unexpected), unless they are depleted. However, for the waitlist group, the effect was absent or partly reversed.

Conclusion

In conclusion, there seems to be group differences at Time 2 and Time 3, between the experimental conditions and the control group. However, the effects in the reflection group appear not only stronger, but also more robust (i.e,. they are the only ones lasting at Time 3). Furthermore, there are also several three-way interactions between implicit attitudes, ego depletion, and group, as expected. The nature of the interactions do not seem however to perfectly align with our original predictions. A deeper exploration of the meaning of these interactions will be required.

Full Code & References

The package references and the full script of executive code contained in this document is reproduced in the tabs below.

…

Package References

sessionInfo() %>% report %>% summary

The analysis was done using the R Statistical language (v4.2.2; R Core Team, 2022) on Windows 10 x64, using the packages iterators (v1.0.14), doParallel (v1.0.17), eefAnalytics (v1.0.6), emmeans (v1.8.5), interactions (v1.1.5), sjlabelled (v1.2.0), performance (v0.10.3), see (v0.7.5), modelbased (v0.8.6), report (v0.5.7.4), foreach (v1.5.2), datawizard (v0.7.1.2), bestNormalize (v1.9.0), missForest (v1.5), rempsyc (v0.1.2.2), visdat (v0.6.0), naniar (v1.0.0), ggplot2 (v3.4.2), dplyr (v1.1.2), tidyr (v1.3.0) and psych (v2.3.3).

report::cite_packages(sessionInfo())

Full Code

library(rempsyc)
library(dplyr)
library(interactions)
library(performance)
library(see)
library(report)
library(datawizard)
library(modelbased)
library(ggplot2)
library(bestNormalize)
library(psych)
library(visdat)
library(missForest)
library(doParallel)
library(ggplot2)
library(emmeans)
library(sjlabelled)
library(eefAnalytics)
library(tidyr)

# Read data
data <- readRDS("Data/finaldataset_n381.rds")

report_participants(data, threshold = 1) %>% cat

# Allocation ratio
report(data$T1_Group)
report(data$T2_Condition)

sessionInfo() %>% report %>% summary

report::cite_packages(sessionInfo())
data %>% 
  filter(is.na(T1_Group)) %>% 
  nrow()

data <- data %>% 
  filter(!is.na(T1_Group))

report_participants(data, threshold = 1) %>% cat

# Allocation ratio
report(data$T1_Group)
report(data$T2_Condition)

data <- data %>% 
  mutate(part.percent = convert_na_to(part.percent, 1))

data %>% 
  filter(part.percent < 2/3) %>% 
  count(T1_Group)

data2 <- data

data <- data %>% 
  filter(part.percent >= 2/3)

report_participants(data, threshold = 1) %>% cat

# Allocation ratio
report(data$T1_Group)
report(data$T2_Condition)

data <- data %>% 
    mutate(att_check = rowSums(
      select(., T1_attention1, T2_attention1, T3_attention1), na.rm = TRUE))

data %>% 
  count(att_check)

data <- data %>% 
  filter(att_check >= 2)

report_participants(data, threshold = 1) %>% cat

report_participants(data, threshold = 1, group = "T1_Group") %>% cat

# Allocation ratio
report(data$T1_Group)
report(data$T2_Condition)

get_label(data$T1_recruitment) %>% cat

report(data$T1_recruitment)

data %>% 
    count(T1_recruitment, sort = TRUE)

data %>% 
  nice_density("age", histogram = TRUE)

data %>% 
    count(gender, sort = TRUE)

get_label(data$T1_psycho.class) %>% cat

data %>% 
    count(T1_psycho.class, sort = TRUE)

get_label(data$T1_virtual.reality) %>% cat

data %>% 
    count(T1_virtual.reality, sort = TRUE)

get_label(data$T1_medi.experience) %>% cat

data %>% 
    count(T1_medi.experience, sort = TRUE, .drop = FALSE)

get_label(data$T1_disorders) %>% cat

data %>% 
    count(T1_disorders, sort = TRUE)

get_label(data$T1_vision) %>% cat

data %>% 
    count(T1_vision, sort = TRUE)

get_label(data$T1_phone) %>% cat

data %>% 
    count(T1_phone, sort = TRUE)

get_label(data$T1_quebec) %>% cat

data %>% 
    count(T1_quebec, sort = TRUE)

get_label(data$T1_student) %>% cat

data %>% 
    count(T1_student, sort = TRUE)

get_label(data$T1_student.program_cat) %>% cat

report(data$T1_student.program)

data %>% 
    count(T1_student.program_cat, sort = TRUE) %>% 
  filter(n > 3)

get_label(data$T1_workplace) %>% cat

report(data$T1_workplace)

data %>% 
    count(T1_workplace_cat, sort = TRUE) %>% 
  filter(n > 1)

get_label(data$T3_post.medipractice) %>% cat

data %>% 
    count(T3_post.medipractice, sort = TRUE)

get_label(data$T3_medipractice.which) %>% cat

report(data$T3_medipractice.which)

data %>% 
    count(T3_medipractice.which, sort = TRUE)

get_label(data$T3_medipractice.time) %>% cat

data %>% 
    count(T3_medipractice.time, sort = TRUE)

get_label(data$T3_choice.medicomp) %>% cat

data %>% 
    count(T3_choice.medicomp, sort = TRUE)

get_label(data$T3_followup) %>% cat

data %>% 
    count(T3_followup, sort = TRUE)

get_label(data$T3_consent) %>% cat

data %>% 
    count(T3_consent, sort = TRUE)

# Check for nice_na
nice_na(data, scales = c(
  "T1_BSCS", "T1_BAQ", "T1_NOBAGS", "T1_attitude", "T1_dehumanization", 
  "T1_WHS", "T1_CLS", "T2_NOBAGS", "T2_attitude", "T2_dehumanization", 
  "T2_SMS5", "T2_PANAS", "T2_WHS", "T2_CLS", "T2_charity", "T3_NOBAGS", 
  "T3_attitude", "T3_dehumanization", "T3_WHS", "T3_CLS"))


# Smaller subset of data for easier inspection
data %>%
  # select(manualworkerId:att_check2_raw, 
  #        condition:condition_dum) %>%
  vis_miss

# Let's use Little's MCAR test to confirm
# We have to proceed by "scale" because the function can only
# support 30 variables max at a time
library(naniar)

# We only check for the variable that had missing data, charity

# Have to divide this one in two because it is too large for the function
data %>% 
  select(T2_charity.moisson1_1:T2_charity.suzuki2_1) %>% 
  mcar_test

data %>% 
  select(T2_charity.conserv1_1:T2_charity.armee2_1) %>% 
  mcar_test

# Need logical and character variables as factors for missForest
# "Error: Can not handle categorical predictors with more than 53 categories."
new.data <- data %>% 
  select(-c(T1_student.program, # T1_student.program = Too many categories (> 53)
            T1_already.participated, # T1_already.participated = lead to error
            T3_consent)) %>% # T3_consent = lead to error
  mutate(across(c(where(is.character), where(is.logical)), as.factor)) %>% 
  as.data.frame()

# Parallel processing
registerDoParallel(cores = 4)

# Variables
set.seed(100)
data.imp <- missForest(new.data, verbose = TRUE, parallelize = "variables")
# Total time is 2 sec (4*0.5) - 4 cores

# Extract imputed dataset
data <- data.imp$ximp

# Reverse code items 2, 4, 6, 7
data <- data %>% 
  mutate(across(contains("BSCS"), .names = "{col}r"))

data <- data %>% 
  mutate(across(ends_with(paste("BSCS", c(2, 4, 6, 7), sep = "_")), 
                ~nice_reverse(.x, 5), .names = "{col}r"))

# Get mean BSCS
data <- data %>% 
  mutate(T1_BSCS = rowMeans(pick(T1_BSCS_1r:T1_BSCS_7r)))

# Get alpha & omega
data %>% 
  select(T1_BSCS_1r:T1_BSCS_7r) %>% 
  omega(nfactors = 1)

# Reverse code item 7
data <- data %>% 
  mutate(across(contains("BAQ"), .names = "{col}r"))

data <- data %>% 
  mutate(across(T1_BAQ_7, ~nice_reverse(.x, 7), .names = "{col}r"))

# Get mean BAQ
data <- data %>% 
  mutate(T1_BAQ = rowMeans(pick(T1_BAQ_1r:T1_BAQ_12r)))

# Get alpha & omega
data %>% 
  select(T1_BAQ_1r:T1_BAQ_12r) %>% 
  omega(nfactors = 1)

data <- data %>% 
  mutate(T1_attitude = rowMeans(pick(T1_attitude_1:T1_attitude_9)),
         T2_attitude = rowMeans(pick(T2_attitude_1:T2_attitude_9)),
         T3_attitude = rowMeans(pick(T3_attitude_1:T3_attitude_9)))

# Get alpha & omega
data %>% 
  select(T1_attitude_1:T1_attitude_9) %>% 
  omega(nfactors = 1)

data %>% 
  select(T2_attitude_1:T2_attitude_9) %>% 
  omega(nfactors = 1)

data %>% 
  select(T3_attitude_1:T3_attitude_9) %>% 
  omega(nfactors = 1)

data <- data %>% 
  mutate(T1_dehumanization = rowMeans(pick(T1_dehumanization_1:T1_dehumanization_9)),
         T2_dehumanization = rowMeans(pick(T2_dehumanization_1:T2_dehumanization_9)),
         T3_dehumanization = rowMeans(pick(T3_dehumanization_1:T3_dehumanization_9)))

# Get alpha & omega
data %>% 
  select(T1_dehumanization_1:T1_dehumanization_9) %>% 
  omega(nfactors = 1)

data %>% 
  select(T2_dehumanization_1:T2_dehumanization_9) %>% 
  omega(nfactors = 1)

data %>% 
  select(T3_dehumanization_1:T3_dehumanization_9) %>% 
  omega(nfactors = 1)

data <- data %>% 
  mutate(across(contains("NOBAGS"), .names = "{col}r"))

# Reverse code NOBAGS (items 1:2, 5:6, 10,12, 14:16, 20)
data <- data %>%
  mutate(across(ends_with(paste0("NOBAGS.", c(
    "1_1", "1_2", "3_1", "3_2", "6_1", "8_1", "10_1", "12_1", "16_1"))), 
    ~nice_reverse(.x, 4), .names = "{col}r"))

# Get mean NOBAGS
data <- data %>% 
  mutate(T1_NOBAGS = rowMeans(pick(T1_NOBAGS.1_1r:T1_NOBAGS.16_1r)),
         T2_NOBAGS = rowMeans(pick(T2_NOBAGS.1_1r:T2_NOBAGS.16_1r)),
         T3_NOBAGS = rowMeans(pick(T3_NOBAGS.1_1r:T3_NOBAGS.16_1r)))

# Get alpha & omega
data %>% 
  select(T1_NOBAGS.1_1r:T1_NOBAGS.16_1r) %>% 
  omega(nfactors = 1)

data %>% 
  select(T2_NOBAGS.1_1r:T2_NOBAGS.16_1r) %>% 
  omega(nfactors = 1)

data %>% 
  select(T3_NOBAGS.1_1r:T3_NOBAGS.16_1r) %>% 
  omega(nfactors = 1)

data <- data %>% 
  mutate(T1_WHS = rowMeans(pick(T1_WHS_1:T1_WHS_6)),
         T2_WHS = rowMeans(pick(T2_WHS_1:T2_WHS_6)),
         T3_WHS = rowMeans(pick(T3_WHS_1:T3_WHS_6)))

# Get alpha & omega
data %>% 
  select(T1_WHS_1:T1_WHS_6) %>% 
  omega(nfactors = 1)

data %>% 
  select(T2_WHS_1:T2_WHS_6) %>% 
  omega(nfactors = 1)

data %>% 
  select(T3_WHS_1:T3_WHS_6) %>% 
  omega(nfactors = 1)

data <- data %>% 
  mutate(T1_CLS = rowMeans(pick(T1_CLS_1:T1_CLS_21)),
         T2_CLS = rowMeans(pick(T2_CLS_1:T2_CLS_21)),
         T3_CLS = rowMeans(pick(T3_CLS_1:T3_CLS_21)))

# Get alpha & omega
data %>% 
  select(T1_CLS_1:T1_CLS_21) %>% 
  omega(nfactors = 1)

data %>% 
  select(T2_CLS_1:T2_CLS_21) %>% 
  omega(nfactors = 1)

data %>% 
  select(T3_CLS_1:T3_CLS_21) %>% 
  omega(nfactors = 1)

data <- data %>% 
  mutate(across(contains("SMS5"), .names = "{col}r"))

# Reverse code SMS5 (items 3 et 5)
data <- data %>%
  mutate(across(ends_with(paste0("SMS5_", c(3, 5))), 
                ~nice_reverse(.x, 5), .names = "{col}r"))

# Get mean SMS5
data <- data %>% 
  mutate(T2_SMS5 = rowMeans(pick(T2_SMS5_1r:T2_SMS5_6r)))

# Get alpha & omega
data %>% 
  select(T2_SMS5_1r:T2_SMS5_6r) %>% 
  omega(nfactors = 1)

# Get mean of PANAS
# Positive affect = 1, 3, 5, 7, 9
# Negative affect = 2, 4, 6, 8, 10
data <- data %>% mutate(
  T2_PANAS_pos = rowMeans(pick(paste0("T2_PANAS_", seq(1, 9, 2)))),
  T2_PANAS_neg = rowMeans(pick(paste0("T2_PANAS_", seq(2, 10, 2)))))

# Get alpha & omega
# Positive affect
data %>% 
  select(all_of(paste0("T2_PANAS_", seq(1, 9, 2)))) %>% 
  omega(nfactors = 1)

# Negative affect
data %>% 
  select(all_of(paste0("T2_PANAS_", seq(2, 10, 2)))) %>% 
  omega(nfactors = 1)

data <- data %>% mutate(
  T2_Charity = rowMeans(pick(contains("charity") & ends_with("1_1"))),
  T2_Familiarity = rowMeans(pick(contains("charity") & ends_with("2_1"))))

# Get alpha & omega
data %>% 
  select(contains("charity") & ends_with("1_1")) %>% 
  omega(nfactors = 1)

# Create new variable blastintensity.duration
data <- data %>% 
  mutate(T1_blastintensity.duration = T1_blastintensity * T1_blastduration,
         T2_blastintensity.duration = T2_blastintensity * T2_blastduration,
         T3_blastintensity.duration = T3_blastintensity * T3_blastduration)

data_temp <- data
data <- data2
data2 <- data
data <- data_temp
# Make list of DVs
col.list <- c("T1_blastintensity", "T1_blastduration", "T1_blastintensity.duration",
              "T2_blastintensity", "T2_blastduration", "T2_blastintensity.duration",
              "T3_blastintensity", "T3_blastduration", "T3_blastintensity.duration",
              "T1_BSCS", "T1_BAQ", "T1_attitude", "T2_attitude", "T3_attitude",
              "T1_dehumanization", "T2_dehumanization", "T3_dehumanization", 
              "T1_NOBAGS", "T1_WHS", "T2_WHS", "T3_WHS", "T1_CLS", "T2_CLS", "T3_CLS", 
              "T2_SMS5", "T2_PANAS_pos", "T2_PANAS_neg", "T2_Charity", 
              "T2_Familiarity", "T2_memory.altruistic", "T2_memory.aggressive")

lapply(col.list, function(x) 
  nice_normality(data, 
                 variable = x, 
                 title = x,
                 group = "T1_Group",
                 shapiro = TRUE,
                 histogram = TRUE))

predict_bestNormalize <- function(var) {
  x <- bestNormalize(var, standardize = FALSE, allow_orderNorm = FALSE)
  print(cur_column())
  print(x$chosen_transform)
  cat("\n")
  y <- predict(x)
  attr(y, "transform") <- c(attributes(y), attributes(x$chosen_transform)$class[1])
  y
}

set.seed(42)
data <- data %>% 
  mutate(across(all_of(col.list), 
                predict_bestNormalize,
                .names = "{.col}.t"))
col.list <- paste0(col.list, ".t")

# Group normality
named.col.list <- setNames(col.list, unlist(lapply(data, function(x) attributes(x)$transform)))
lapply(named.col.list, function(x) 
  nice_normality(data, 
                 x, 
                 "T1_Group",
                 shapiro = TRUE,
                 title = x,
                 histogram = TRUE))

# Plotting variance
plots(lapply(col.list, function(x) {
  nice_varplot(data, x, group = "T1_Group")
  }),
  n_columns = 2)

plots(lapply(col.list, function(x) {
  plot_outliers(data, x, group = "T1_Group", ytitle = x, binwidth = 0.3)
  }),
  n_columns = 2)

data %>% 
  as.data.frame %>% 
  filter(T1_Group == "Waitlist") %>% 
  find_mad(col.list, criteria = 3)

data %>% 
  as.data.frame %>% 
  filter(T1_Group == "Reflection") %>% 
  find_mad(col.list, criteria = 3)

data %>% 
  as.data.frame %>% 
  filter(T1_Group == "Meditation") %>% 
  find_mad(col.list, criteria = 3)

data.na <- na.omit(data[col.list])
x <- check_outliers(data.na[-c(15:16)], method = "mcd", threshold = 500)
# We have to exclude two variables that are too large following the transformations
# Otherwise we get an error:
# Error in solve.default(cov, ...) : 
#  system is computationally singular: reciprocal condition number = 8.99059e-20
x

plot(x)

# Winsorize variables of interest with MAD
data <- data %>% 
  group_by(T1_Group) %>% 
  mutate(across(all_of(col.list), 
                winsorize_mad,
                .names = "{.col}.w")) %>% 
  ungroup()

# Update col.list
col.list <- paste0(col.list, ".w")

data <- data %>%
  mutate(across(all_of(col.list), standardize, .names = "{col}.s"))

# Update col.list
col.list <- paste0(col.list, ".s")

# Let's replace original variables with the transformed variables
data[gsub(".t.w.s", "", col.list)] <- data[col.list]

# If we decide to only center variables instead of standardizing them
# (as in the preregistration), then let's do this instead
data <- data %>%
  mutate(across(all_of(gsub(".t.w.s", "", col.list)), 
                \(x) standardize(x, center = TRUE, scale = FALSE), 
                .names = "{col}"))

data_temp <- data
data <- data2
col.list_temp <- col.list
col.list <- c("T1_blastintensity", "T1_blastduration", "T1_blastintensity.duration",
              "T2_blastintensity", "T2_blastduration", "T2_blastintensity.duration",
              "T3_blastintensity", "T3_blastduration", "T3_blastintensity.duration",
              "T1_BSCS", "T1_BAQ", "T1_attitude", "T2_attitude", "T3_attitude",
              "T1_dehumanization", "T2_dehumanization", "T3_dehumanization", 
              "T1_NOBAGS", "T1_WHS", "T2_WHS", "T3_WHS", "T1_CLS", "T2_CLS", "T3_CLS", 
              "T2_SMS5", "T2_PANAS_pos", "T2_PANAS_neg", "T2_Charity",
              "T2_memory.altruistic", "T2_memory.aggressive")
data2 <- data
data <- data_temp
col.list <- col.list_temp 
# Specify the order of factor levels for "Group". 
# Otherwise R will alphabetize them.
data$T1_Group <- factor(data$T1_Group, levels = c("Meditation", "Reflection", "Waitlist"))

# Define our dependent variables
DV <- data %>% select(T2_NOBAGS:T2_Charity) %>% names

# First column (which variable)
Variable <- rep(DV, each = 3)

# Second column (which comparison)
Comparison <- rep(c("MeditationvsCTR", 
                    "ReflectionvsCTR", 
                    "MeditationvsReflection"), 
                  length(DV))
# 14 == number of DV

# Make list of all formulas
formulas <- c(
  "T2_NOBAGS ~ T1_Group * T1_NOBAGS",
  "T2_attitude ~ T1_Group * T1_attitude",
  "T2_dehumanization ~ T1_Group * T1_dehumanization",
  "T2_IAT ~ T1_Group * T1_IAT",
  "T2_SMS5 ~ T1_Group",
  "T2_PANAS_pos ~ T1_Group",
  "T2_PANAS_neg ~ T1_Group",
  "T2_blastintensity ~ T1_Group * T1_blastintensity",
  "T2_blastduration ~ T1_Group * T1_blastduration",
  "T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration",
  "T2_memory.altruistic ~ T1_Group",
  "T2_memory.aggressive ~ T1_Group",
  "T2_WHS ~ T1_Group * T1_WHS",
  "T2_CLS ~ T1_Group * T1_CLS",
  "T2_Charity ~ T1_Group * T2_Familiarity" 
)

# Make list of all models
models.list <- sapply(formulas, lm, data = data, simplify = FALSE, USE.NAMES = TRUE)

## Attempt with nice_lm_contrasts 
x <- nice_lm_contrasts(models.list, group = "T1_Group", data = data)

x.table <- nice_table(x, highlight = TRUE)
x.table
# Make list of all formulas
formulas2 <- c(
  "T3_NOBAGS ~ T1_Group * T1_NOBAGS",
  "T3_attitude ~ T1_Group * T1_attitude",
  "T3_dehumanization ~ T1_Group * T1_dehumanization",
  "T3_IAT ~ T1_Group * T1_IAT",
  "T3_blastintensity ~ T1_Group * T1_blastintensity",
  "T3_blastduration ~ T1_Group * T1_blastduration",
  "T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration",
  "T3_WHS ~ T1_Group * T1_WHS",
  "T3_CLS ~ T1_Group * T1_CLS"
)

# Make list of all models
models.list2 <- sapply(formulas2, lm, data = data, simplify = FALSE, USE.NAMES = TRUE)

## Attempt with nice_lm_contrasts 
x2 <- nice_lm_contrasts(models.list2, group = "T1_Group", data = data)

x.table2 <- nice_table(x2, highlight = TRUE)
x.table2
nice_violin(data, 
            group = "T1_Group", 
            response = "T2_attitude",
            obs = TRUE,
            signif_annotation = c("*"),
            signif_yposition = 1.5,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(models.list$`T2_attitude ~ T1_Group * T1_attitude`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Attitude", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(models.list$`T2_attitude ~ T1_Group * T1_attitude`)
plot(contrasts, estimate_means(models.list$`T2_attitude ~ T1_Group * T1_attitude`)) +
  ylab(paste("Adjusted", "Attitude", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_PANAS_pos",
            obs = TRUE,
            comp1 = 1,
            comp2 = 3,
            has.d = TRUE,
            d.x = 1.75,
            d.y = 2)

means <- estimate_means(models.list$`T2_PANAS_pos ~ T1_Group`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Positive Affect", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(models.list$`T2_PANAS_pos ~ T1_Group`)
plot(contrasts, estimate_means(models.list$`T2_PANAS_pos ~ T1_Group`)) +
  ylab(paste("Adjusted", "Positive Affect", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_blastintensity",
            obs = TRUE,
            signif_annotation = c("**", "***"),
            signif_yposition = c(3.5, 4),
            signif_xmin = c(1, 2),
            signif_xmax = c(2, 3))

means <- estimate_means(models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)
plot(contrasts, estimate_means(
  models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)) +
  ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_blastduration",
            obs = TRUE,
            signif_annotation = c("**", "***"),
            signif_yposition = c(3.5, 4),
            signif_xmin = c(1, 2),
            signif_xmax = c(2, 3))

means <- estimate_means(models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Duration", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)
plot(contrasts, estimate_means(
  models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)) +
  ylab(paste("Adjusted", "Blast Duration", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_blastintensity.duration",
            obs = TRUE,
            signif_annotation = c("*", "**"),
            signif_yposition = c(3.5, 4),
            signif_xmin = c(1, 2),
            signif_xmax = c(2, 3))

means <- estimate_means(
  models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
plot(contrasts, estimate_means(
  models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_memory.altruistic",
            obs = TRUE,
            signif_annotation = c("*"),
            signif_yposition = 3.5,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(
  models.list$`T2_memory.altruistic ~ T1_Group`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "ms to remember altruistic event", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list$`T2_memory.altruistic ~ T1_Group`)
plot(contrasts, estimate_means(
  models.list$`T2_memory.altruistic ~ T1_Group`)) +
  ylab(paste("Adjusted", "ms to remember altruistic event", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T2_CLS",
            obs = TRUE,
            signif_annotation = c("*", "***"),
            signif_yposition = c(3.2, 2.5),
            signif_xmin = c(1, 2),
            signif_xmax = c(3, 3))

means <- estimate_means(models.list$`T2_CLS ~ T1_Group * T1_CLS`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Compassion", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(models.list$`T2_CLS ~ T1_Group * T1_CLS`)
plot(contrasts, estimate_means(models.list$`T2_CLS ~ T1_Group * T1_CLS`)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_blastintensity",
            obs = TRUE,
            signif_annotation = c("**"),
            signif_yposition = 2.5,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)
plot(contrasts, estimate_means(
  models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)) +
  ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_blastduration",
            obs = TRUE,
            signif_annotation = c("**"),
            signif_yposition = 3,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Duration", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)
plot(contrasts, estimate_means(
  models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)) +
  ylab(paste("Adjusted", "Blast Duration", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_blastintensity.duration",
            obs = TRUE,
            signif_annotation = c("*"),
            signif_yposition = 3,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(
  models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
plot(contrasts, estimate_means(
  models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_WHS",
            obs = TRUE,
            signif_annotation = c("**"),
            signif_yposition = 2.5,
            signif_xmin = 2,
            signif_xmax = 3)

means <- estimate_means(models.list2$`T3_WHS ~ T1_Group * T1_WHS`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Willingness to Help", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(
  models.list2$`T3_WHS ~ T1_Group * T1_WHS`)
plot(contrasts, estimate_means(
  models.list2$`T3_WHS ~ T1_Group * T1_WHS`)) +
  ylab(paste("Adjusted", "Willingness to Help", "Mean")) +
  theme_modern()

nice_violin(data, 
            group = "T1_Group", 
            response = "T3_CLS",
            obs = TRUE,
            signif_annotation = c("*", "**"),
            signif_yposition = c(3.2, 2.7),
            signif_xmin = c(1, 2),
            signif_xmax = c(3, 3))

means <- estimate_means(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)

ggplot(means, aes(x = T1_Group, y = Mean)) +
  geom_line(aes(group = 1)) +
  geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
  ylab(paste("Adjusted", "Compassion", "Mean")) +
  theme_modern()

contrasts <- estimate_contrasts(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)
plot(contrasts, estimate_means(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)) +
  ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
  theme_modern()

data3 <- data2 %>% 
  mutate(part.percent = ifelse(T1_Group == "Waitlist", 0, part.percent),
         part.percent = part.percent * 100,
         T2_memory.altruistic = T2_memory.altruistic * -1,
         T1_GroupReflection = as.numeric(T1_Group == "Reflection"),
         T1_GroupMeditation = as.numeric(T1_Group == "Meditation"),
         across(contains("blast"), \(x) x * -1)) %>% 
  as.data.frame()

caceOutput <- caceSRTBoot(
  T2_attitude ~ T1_attitude + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Attitude (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T2_PANAS_pos ~ T1_GroupMeditation,
  intervention = "T1_GroupMeditation",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Positive Affect (Meditation VS Waitlist)")

caceOutput <- caceSRTBoot(
  T2_blastintensity ~ T1_blastintensity + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Intensity (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T2_blastduration ~ T1_blastduration + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Duration (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T2_blastintensity.duration ~ T1_blastintensity.duration + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Intensity * Duration (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T2_memory.altruistic ~ T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "ms to remember altruistic event (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T2_CLS ~ T1_CLS + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T2_CLS ~ T1_CLS + T1_GroupMeditation,
  intervention = "T1_GroupMeditation",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Meditation VS Waitlist)")

caceOutput <- caceSRTBoot(
  T3_attitude ~ T1_attitude + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Attitude (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T3_blastintensity ~ T1_blastintensity + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Intensity (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T3_blastduration ~ T1_blastduration + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Duration (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T3_blastintensity.duration ~ T1_blastintensity.duration + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Blast Intensity * Duration (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T3_WHS ~ T1_WHS + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Willingness to Help (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T3_CLS ~ T1_CLS + T1_GroupReflection,
  intervention = "T1_GroupReflection",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Reflection VS Waitlist)")

caceOutput <- caceSRTBoot(
  T3_CLS ~ T1_CLS + T1_GroupMeditation,
  intervention = "T1_GroupMeditation",
  compliance = "part.percent",
  nBoot = 2000,
  data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Meditation VS Waitlist)")

# CREATE OUR DUMMY VARIABLES FOR T1_Group!

data$T1_GroupReflection <- as.numeric(data$T1_Group == "Reflection")
data$T1_GroupMeditation <- as.numeric(data$T1_Group == "Meditation")

################################################
################################################

# T2_blastintensity.duration - IAT
T2_blastintensity.duration.IAT <- lm(
  T2_blastintensity.duration ~ T1_GroupReflection + T1_GroupMeditation +
    T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + T2_IAT:T2_Condition + 
    T1_GroupReflection:T2_IAT:T2_Condition + T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_blastintensity.duration.IAT)

T2_blastintensity.duration.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# Make interaction plot
interact_plot(T2_blastintensity.duration.IAT, pred = "T2_IAT", modx = "T2_Condition",
              interval = TRUE, mod2 = "T1_GroupMeditation",
              x.label = "Implicit Aggression (IAT)",
              y.label = "Blast Intensity * Duration (Taylor)", 
              mod2.labels = c("Waitlist Group", "Meditation Group"),
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

# T2_blastintensity - IAT
T2_blastintensity.IAT <- lm(
  T2_blastintensity ~ T1_GroupReflection + T1_GroupMeditation +
    T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + 
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_blastintensity.IAT)

T2_blastintensity.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# T2_blastduration - IAT
T2_blastduration.IAT <- lm(
  T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
    T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_blastduration.IAT)

T2_blastduration.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# Make interaction plot
interact_plot(T2_blastduration.IAT, pred = "T2_IAT", modx = "T2_Condition",
              interval = TRUE, mod2 = "T1_GroupMeditation",
              x.label = "Implicit Aggression (IAT)",
              y.label = "Blast Duration (Taylor)", 
              mod2.labels = c("Waitlist Group", "Meditation Group"),
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

# T2_Charity - IAT
T2_Charity.IAT <- lm(
  T2_Charity ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
    T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +  
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition + T2_Familiarity, 
  data = data)

check_model(T2_Charity.IAT)

T2_Charity.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# Compassion - IAT
T2_CLS.IAT <- lm(
  T2_CLS ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + 
    T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_CLS.IAT)

T2_CLS.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# Make interaction plot
interact_plot(T2_CLS.IAT, pred = "T2_IAT", modx = "T2_Condition",
              interval = TRUE, mod2 = "T1_GroupMeditation",
              x.label = "Implicit Aggression (IAT)",
              y.label = "Compassionate Love",
              mod2.labels = c("Waitlist Group", "Meditation Group"),
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

# T2_WHS - IAT
T2_WHS.IAT <- lm(
  T2_WHS ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + T2_Condition + 
    T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT + T1_GroupReflection:T2_Condition +
    T1_GroupMeditation:T2_Condition + T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_WHS.IAT)

T2_WHS.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# T2_memory.altruistic - IAT
T2_memory.altruistic.IAT <- lm(
  T2_memory.altruistic ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + 
    T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +  
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_memory.altruistic.IAT)

T2_memory.altruistic.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# Make interaction plot
interact_plot(T2_memory.altruistic.IAT, pred = "T2_IAT", modx = "T2_Condition",
              interval = TRUE, mod2 = "T1_GroupMeditation",
              x.label = "Implicit Aggression (IAT)",
              y.label = "Reaction Time (Altruistic Memory)",
              mod2.labels = c("Waitlist Group", "Meditation Group"),
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

# T2_memory.aggressive - IAT
T2_memory.aggressive.IAT <- lm(
  T2_memory.aggressive ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + 
    T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +  
    T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition + 
    T1_GroupMeditation:T2_IAT:T2_Condition, data = data)

check_model(T2_memory.aggressive.IAT)

T2_memory.aggressive.IAT %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# T2_blastintensity.duration - NOBAGS
T2_blastintensity.duration.NOBAGS <- lm(
  T2_blastintensity.duration ~ T1_GroupReflection + T1_GroupMeditation +
    T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + 
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_blastintensity.duration.NOBAGS)

T2_blastintensity.duration.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# T2_blastintensity - NOBAGS
T2_blastintensity.NOBAGS <- lm(
  T2_blastintensity ~ T1_GroupReflection + T1_GroupMeditation +
    T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + 
    T1_GroupMeditation:T2_NOBAGS +  T1_GroupReflection:T2_Condition + 
    T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition + 
    T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_blastintensity.NOBAGS)

T2_blastintensity.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# T2_blastduration - NOBAGS
T2_blastduration.NOBAGS <- lm(
  T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
    T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + 
    T1_GroupMeditation:T2_NOBAGS +  T1_GroupReflection:T2_Condition +
    T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition +
    T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_blastduration.NOBAGS)

T2_blastduration.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)

# Make interaction plot
interact_plot(T2_blastduration.NOBAGS, pred = "T2_NOBAGS", modx = "T2_Condition",
              interval = TRUE, 
              x.label = "Normative beliefs about aggression (NOBAGS)",
              y.label = "Blast Duration (Taylor)",
              legend.main = "Condition") +
  theme(text = element_text(size = 25))

# T2_blastduration.NOBAGS <- lm(
#   T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
#     T2_NOBAGS + T2_Condition_dum + T1_GroupReflection:T2_NOBAGS + 
#     T1_GroupMeditation:T2_NOBAGS +  T1_GroupReflection:T2_Condition_dum +
#     T1_GroupMeditation:T2_Condition_dum + T2_NOBAGS:T2_Condition_dum +
#     T1_GroupReflection:T2_NOBAGS:T2_Condition_dum + 
#     T1_GroupMeditation:T2_NOBAGS:T2_Condition_dum, data = data)

# T2_blastduration.NOBAGS %>%
#   nice_lm_slopes(predictor = "T2_Condition_dum",
#                  moderator = "T2_NOBAGS") %>% 
#   nice_table(highlight = TRUE)
# 
# ?nice_slopes
# 
# probe_interaction(T2_blastintensity.duration.NOBAGS, pred = T2_NOBAGS, 
#                   modx = T2_Condition, digits = 3)
# 
# emtrends(T2_blastintensity.duration.NOBAGS, pairwise ~ T2_Condition, var = "T2_NOBAGS")

# T2_Charity - NOBAGS
T2_Charity.NOBAGS <- lm(
  T2_Charity ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS +
    T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition + T2_Familiarity, 
  data = data)

check_model(T2_Charity.NOBAGS)

T2_Charity.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# Compassion - NOBAGS
T2_CLS.NOBAGS <- lm(
  T2_CLS ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS + 
    T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +  
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_CLS.NOBAGS)

T2_CLS.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# T2_WHS - NOBAGS
T2_WHS.NOBAGS <- lm(
  T2_WHS ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS + T2_Condition + 
    T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS + 
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + 
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_WHS.NOBAGS)

T2_WHS.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# T2_memory.altruistic - NOBAGS
T2_memory.altruistic.NOBAGS <- lm(
  T2_memory.altruistic ~ T1_GroupReflection + T1_GroupMeditation +
    T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + 
    T1_GroupMeditation:T2_NOBAGS + T1_GroupReflection:T2_Condition + 
    T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition + 
    T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_memory.altruistic.NOBAGS)

T2_memory.altruistic.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
# T2_memory.aggressive - NOBAGS
T2_memory.aggressive.NOBAGS <- lm(
  T2_memory.aggressive ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS + 
    T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +  
    T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + 
    T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition + 
    T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)

check_model(T2_memory.aggressive.NOBAGS)

T2_memory.aggressive.NOBAGS %>% 
  nice_lm() %>% 
  nice_table(highlight = TRUE)
data %>%
  select(starts_with("T1_attitude")) %>% 
  get_label %>% 
  lapply(function(x) gsub(".*- ", "", x)) %>% 
  unlist() %>% unname

social.groups <- c("Blacks", "Homeless", "Native", "Muslims", "Refugees", "Women",
                   "Animals", "Elderly", "Whites")

charities <- data %>%
  select(ends_with("1_1") & contains("charity")) %>% 
  get_label %>% 
  lapply(function(x) gsub(".*- ", "", x)) %>% 
  unlist() %>% unname
charities

regions <- c("Montreal", "Quebec", "Canada", "International")
regions

data %>%
  select(T1_attitude_1:T1_attitude_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Positive Explicit Attitude",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

data %>%
  select(T1_dehumanization_1:T1_dehumanization_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Humanization",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

data %>%
  select(T2_attitude_1:T2_attitude_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Positive Explicit Attitude",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

data %>%
  select(T2_dehumanization_1:T2_dehumanization_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Humanization",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

data %>%
  select(contains("charity") & ends_with("1_1")) %>% 
  pivot_longer(cols = everything(), 
               names_to = "charity",
               values_to = "donation") %>% 
  mutate(charity = factor(charity, labels = charities)) %>% 
  nice_violin(group = "charity",
              response = "donation",
              ytitle = "Amount Donated",
              CIcap.width = 0.5,
              obs = "jitter",
              order.groups = "increasing",
              xlabels.angle = 75)

data %>%
  select(contains("charity") & ends_with("2_1")) %>% 
  pivot_longer(cols = everything(), 
               names_to = "charity",
               values_to = "familiarity") %>% 
  mutate(charity = factor(charity, labels = charities)) %>% 
  nice_violin(group = "charity",
              response = "familiarity",
              ytitle = "Familiarity with Charity",
              CIcap.width = 0.5,
              obs = "jitter",
              order.groups = "increasing",
              xlabels.angle = 75)

data %>%
  nice_scatter(predictor = "T2_Familiarity",
               response = "T2_Charity",
               ytitle = "Donation Amount",
               xtitle = "Familiarity with Charity",
               has.jitter = TRUE,
               has.legend = TRUE,
               has.r = TRUE,
               has.p = TRUE)

data %>%
  select(contains("charity") & ends_with("1_1")) %>% 
  pivot_longer(cols = everything(), 
               names_to = "charity",
               values_to = "donation") %>% 
  mutate(charity = factor(charity, labels = charities),
         region = case_match(
           charity,
           charities[1:6] ~ regions[1],
           charities[7:12] ~ regions[2],
           charities[13:18] ~ regions[3],
           charities[19:24] ~ regions[4]
         )) %>% 
  nice_violin(group = "region",
              response = "donation",
              ytitle = "Amount Donated",
              obs = "jitter",
              order.groups = "increasing")

data %>%
  select(contains("charity") & ends_with("2_1")) %>% 
  pivot_longer(cols = everything(), 
               names_to = "charity",
               values_to = "familiarity") %>% 
  mutate(charity = factor(charity, labels = charities),
         region = case_match(
           charity,
           charities[1:6] ~ regions[1],
           charities[7:12] ~ regions[2],
           charities[13:18] ~ regions[3],
           charities[19:24] ~ regions[4]
         )) %>% 
  nice_violin(group = "region",
              response = "familiarity",
              ytitle = "Familiarity with Charity",
              obs = "jitter",
              order.groups = "increasing")

data %>%
  select(T3_attitude_1:T3_attitude_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Positive Explicit Attitude",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")

data %>%
  select(T3_dehumanization_1:T3_dehumanization_9) %>% 
  pivot_longer(cols = everything(), 
               names_to = "Group",
               values_to = "attitude") %>% 
  mutate(Group = factor(Group, labels = social.groups)) %>% 
  nice_violin(group = "Group",
              response = "attitude",
              ytitle = "Humanization",
              CIcap.width = 0.3,
              obs = "jitter",
              order.groups = "increasing")
LS0tDQp0aXRsZTogJyoqTG92aW5nLUtpbmRuZXNzIE1lZGl0YXRpb24sIEF0dGl0dWRlcywgUHJvc29jaWFsIEJlaGF2aW91ciwgYW5kIEVnbyBEZXBsZXRpb24qKicNCnN1YnRpdGxlOiAiSXMgdGhlIE1pbmQgTW9yZSBQb3dlcmZ1bCBUaGFuIHRoZSBIZWFydD8gQSBMb29rIGF0IFR3byBMb3ZpbmctS2luZG5lc3MgSW50ZXJ2ZW50aW9ucyINCmF1dGhvcjogIlLDqW1pIFRow6lyaWF1bHQiDQpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCINCm91dHB1dDoNCiAgcm1hcmtkb3duOjpodG1sX2RvY3VtZW50Og0KICAgICAgICB0aGVtZTogY2VydWxlYW4NCiAgICAgICAgaGlnaGxpZ2h0OiBweWdtZW50cw0KICAgICAgICB0b2M6IHllcw0KICAgICAgICB0b2NfZGVwdGg6IDINCiAgICAgICAgdG9jX2Zsb2F0OiB5ZXMNCiAgICAgICAgbnVtYmVyX3NlY3Rpb25zOiBubw0KICAgICAgICBkZl9wcmludDoga2FibGUNCiAgICAgICAgY29kZV9mb2xkaW5nOiBzaG93ICMgb3I6IGhpZGUNCiAgICAgICAgY29kZV9kb3dubG9hZDogeWVzDQogICAgICAgIGFuY2hvcl9zZWN0aW9uczoNCiAgICAgICAgICBzdHlsZTogc3ltYm9sDQogICAgICANCi0tLQ0KDQpgYGB7ciBzZXR1cCwgd2FybmluZz1GQUxTRSwgbWVzc2FnZT1UUlVFLCBlY2hvPUZBTFNFfQ0KZmFzdCA8LSBGQUxTRSAjIE1ha2UgdGhpcyB0cnVlIHRvIHNraXAgdGhlIHRpbWUtY29uc3VtaW5nIGNodW5rcw0KDQp0cmFuc2Zvcm0gPC0gRkFMU0UgIyBNYWtlIHRoaXMgRkFMU0UgdG8gc2tpcCB0aGUgdHJhbnNmb3JtYXRpb24gY2h1bmtzDQoNCmtuaXRyOjpvcHRzX2NodW5rJHNldChkcGkgPSAxMDApDQpgYGANCg0KYGBge3Iga2xpcHB5LCBlY2hvPUZBTFNFLCBpbmNsdWRlPVRSVUV9DQprbGlwcHk6OmtsaXBweShwb3NpdGlvbiA9IGMoJ3RvcCcsICdyaWdodCcpKQ0KYGBgDQoNCiMgSW50cm9kdWN0aW9uDQoNClRoaXMgcmVwb3J0IGRlc2NyaWJlcyB0aGUgcmVzdWx0cyBvZiBhIHByZXJlZ2lzdGVyZWQgc3R1ZHkgYXZhaWxhYmxlIGF0OiBodHRwczovL29zZi5pby9na2Q4cy4NCg0KLS0tDQpOb3RlIGFsc28gdGhhdCB0aGlzIGRhdGEgaGFzIGJlZW4gY2xlYW5lZCBiZWZvcmVoYW5kLiBTZXZlcmFsIGRhdGFzZXRzIG92ZXIgdGhyZWUgbWVhc3VyZW1lbnQgdGltZXMgd2VyZSBtZXJnZWQgKGpvaW5lZCkgdGhyb3VnaCBhbiBpbm5lciBqb2luIHNvIGFzIHRvIGtlZXAgb25seSBwYXJ0aWNpcGFudHMgd2hvIGF0IGxlYXN0IHBhcnRpY2lwYXRlZCBhdCBlYWNoIHN0ZXAgb2YgdGhlIHN0dWR5LiBNaXNzaW5nIGRhdGEgd2lsbCBiZSBpbXB1dGVkIGxhdGVyIG9uLiBEdXBsaWNhdGVzIHdlcmUgYWRkcmVzc2VkIHdpdGggdGhlIGByZW1wc3ljOjpiZXN0X2R1cGxpY2F0ZWAgZnVuY3Rpb24sIHdoaWNoIGtlZXBzIHRoZSBkdXBsaWNhdGUgd2l0aCB0aGUgbGVhc3QgYW1vdW50IG9mIG1pc3NpbmcgdmFsdWVzLCBhbmQgaW4gY2FzZSBvZiB0aWVzLCB0YWtlcyB0aGUgZmlyc3Qgb2NjdXJyZW5jZS4gSG93ZXZlciwgZm9yIGR1cGxpY2F0ZSBwYXJ0aWNpcGF0aW9uIGluIHRoZSBhY3Rpdml0aWVzIGFuZCBleGVyY2lzZXMsIHJhdGhlciB0aGFuIHRoZSBmaXJzdCBvY2N1cnJlbmNlLCB0aGUgb2NjdXJyZW5jZSB3aXRoIHRoZSBoaWdoZXIgY29tcGxldGlvbiBwZXJjZW50YWdlICglIG9mIHRvdGFsIGFjdGl2aXR5IHRpbWUpIHdhcyB0YWtlbiBpbnN0ZWFkLg0KDQojIFBhY2thZ2VzICYgRGF0YSB7LnRhYnNldH0NCg0KIyMgLi4uDQoNCiMjIFBhY2thZ2VzDQoNCmBgYHtyIHBhY2thZ2VzLCB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCByZXN1bHRzPSdhc2lzJ30NCmxpYnJhcnkocmVtcHN5YykNCmxpYnJhcnkoZHBseXIpDQpsaWJyYXJ5KGludGVyYWN0aW9ucykNCmxpYnJhcnkocGVyZm9ybWFuY2UpDQpsaWJyYXJ5KHNlZSkNCmxpYnJhcnkocmVwb3J0KQ0KbGlicmFyeShkYXRhd2l6YXJkKQ0KbGlicmFyeShtb2RlbGJhc2VkKQ0KbGlicmFyeShnZ3Bsb3QyKQ0KbGlicmFyeShiZXN0Tm9ybWFsaXplKQ0KbGlicmFyeShwc3ljaCkNCmxpYnJhcnkodmlzZGF0KQ0KbGlicmFyeShtaXNzRm9yZXN0KQ0KbGlicmFyeShkb1BhcmFsbGVsKQ0KbGlicmFyeShnZ3Bsb3QyKQ0KbGlicmFyeShlbW1lYW5zKQ0KbGlicmFyeShzamxhYmVsbGVkKQ0KbGlicmFyeShlZWZBbmFseXRpY3MpDQpsaWJyYXJ5KHRpZHlyKQ0KDQpgYGANCg0KIyMgRGF0YQ0KDQpgYGB7ciBkYXRhLCB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPVRSVUUsIHJlc3VsdHM9J2FzaXMnfQ0KIyBSZWFkIGRhdGENCmRhdGEgPC0gcmVhZFJEUygiRGF0YS9maW5hbGRhdGFzZXRfbjM4MS5yZHMiKQ0KDQpyZXBvcnRfcGFydGljaXBhbnRzKGRhdGEsIHRocmVzaG9sZCA9IDEpICU+JSBjYXQNCg0KIyBBbGxvY2F0aW9uIHJhdGlvDQpyZXBvcnQoZGF0YSRUMV9Hcm91cCkNCnJlcG9ydChkYXRhJFQyX0NvbmRpdGlvbikNCg0KYGBgDQoNCiMgRGF0YSBDbGVhbmluZw0KDQpgYGB7ciBEYXRhIENsZWFuaW5nLCBjaGlsZD1pZiAoZmFzdCA9PSBGQUxTRSkgJzBfY2xlYW5pbmcuUm1kJywgZXZhbCA9IFRSVUV9DQpgYGANCg0KIyBBc3N1bXB0aW9ucyB7LnRhYnNldH0NCg0KYGBge3IgQXNzdW1wdGlvbnMsIGNoaWxkPWlmIChmYXN0ID09IEZBTFNFKSAnMV9hc3N1bXB0aW9ucy5SbWQnLCBldmFsID0gVFJVRX0NCmBgYA0KDQojIENvbnRyYXN0cyAoR3JvdXAgRGlmZmVyZW5jZXMpIHsudGFic2V0fQ0KDQpgYGB7ciBDb250cmFzdHMsIGNoaWxkPWlmIChmYXN0ID09IEZBTFNFKSAnMl9jb250cmFzdHMuUm1kJywgZXZhbCA9IFRSVUV9DQpgYGANCg0KIyBNb2RlcmF0aW9ucyAoRWdvIERlcGxldGlvbikgey50YWJzZXR9DQoNCiMjIC4uLg0KDQpgYGB7ciBNb2RlcmF0aW9ucywgY2hpbGQ9aWYgKGZhc3QgPT0gRkFMU0UpICczX21vZGVyYXRpb25zLlJtZCcsIGV2YWwgPSBUUlVFfQ0KYGBgDQoNCiMgQXR0aXR1ZGVzICYgQ2hhcml0eSB7LnRhYnNldH0NCg0KIyMgLi4uDQoNCmBgYHtyIENoYXJpdHksIGNoaWxkPWlmIChmYXN0ID09IEZBTFNFKSAnNF9jaGFyaXR5LlJtZCcsIGV2YWwgPSBUUlVFfQ0KYGBgDQoNCiMgRGlzY3Vzc2lvbg0KDQpJbiB0aGlzIHJlcG9ydCwgd2UgYWltZWQgdG8gY29tcGFyZSB0d28gdHlwZXMgb2YgbG92aW5nLWtpbmRuZXNzIG1lZGl0YXRpb24sIG9uZSBtb3JlIGVtYm9kaWVkLCBiYXNlZCBvbiBtZWRpdGF0aW9uLCBhbmQgb25lIG1vcmUgY29nbml0aXZlLCBiYXNlZCBvbiBpbnRlbGxlY3R1YWwgcmVmbGVjdGlvbiwgdG8gYSB3YWl0bGlzdCBjb250cm9sIGdyb3VwLiBXZSBjb21wYXJlZCB0aG9zZSBncm91cHMgb24gc2V2ZXJhbCB2YXJpYWJsZXMgcmVsYXRpbmcgdG8gcHJvc29jaWFsaXR5IChpLmUuLCBvbiBhZmZlY3QsIGF0dGl0dWRlIGFuZCBiZWhhdmlvdXIpLiBHcm91cHMgd2VyZSBtZWFzdXJlZCB0aHJlZSB0aW1lczogV2VlayAwIChUMSksIFdlZWsgNiAoVDIpLCBhbmQgV2VlayAxMyAoVDMpIHNvIGl0IHdhcyBhYmxlIHRvIGNvbXBhcmUgZm9yIGJhc2VsaW5lIGJ1dCBhbHNvIHNlZSBob3cgcm9idXN0IHRoZSBlZmZlY3RzLCBpZiBhbnksIGFyZSB0aHJvdWdoIHRpbWVzLiBXZSB3ZXJlIGFsc28gaW50ZXJlc3RlZCBpbiBhc3Nlc3Npbmcgd2hldGhlciB0aGVzZSBlZmZlY3RzIGRlcGVuZCBvbiBvdGhlciBwZXJzb25hbGl0eSB2YXJpYWJsZXMgKGkuZS4sIG1vZGVyYXRvcnMpLg0KDQojIyMjIEdyb3VwIERpZmZlcmVuY2VzIGF0IFRpbWUgMg0KDQpPdXIgY29udHJhc3RzIGFuYWx5c2VzIGZpcnN0IHJldmVhbGVkIGdyb3VwIGRpZmZlcmVuY2VzIGF0IFRpbWUgMi4gQm90aCB0aGUgbWVkaXRhdGlvbiBhbmQgcmVmbGVjdGlvbiBncm91cHMgc2hvd2VkIG1vZGVyYXRlbHkgbW9yZSBjb21wYXNzaW9uYXRlIGxvdmUgdGhhbiB0aGUgd2FpdGxpc3QgZ3JvdXAsIGJ1dCBvbmx5IHRoZSBtZWRpdGF0aW9uIGdyb3VwIHNob3dlZCBtb2RlcmF0ZWx5IG1vcmUgcG9zaXRpdmUgYWZmZWN0IHRoYW4gdGhlIHdhaXRsaXN0IGdyb3VwLiBIb3dldmVyLCB0aGUgcmVmbGVjdGlvbiBncm91cCBzaG93ZWQgYSBsaXR0bGUgbW9yZSBwb3NpdGl2ZSBleHBsaWNpdCBhdHRpdHVkZXMgdG93YXJkIHZhcmlvdXMgc29jaWFsIGdyb3VwcywgYXMgd2VsbCBhcyBtb2RlcmF0ZWx5IHNob3J0ZXIgcmVhY3Rpb24gdGltZXMgdG8gcmVtZW1iZXIgYW4gYWx0cnVpc3RpYyBldmVudCB0aGFuIHRoZSB3YWl0bGlzdCBncm91cCAoc3VnZ2VzdGluZyB0aGF0IGFsdHJ1aXNtIHdhcyBtb3JlIGNvZ25pdGl2ZWx5IGFjY2Vzc2libGUgdG8gdGhlbSkuIEZ1cnRoZXJtb3JlLCB0aGUgcmVmbGVjdGlvbiBncm91cCBzaG93ZWQgYSBsaXR0bGUgbG93ZXIgYmVoYXZpb3VyYWwgYWdncmVzc2lvbiAoYmxhc3QgaW50ZW5zaXR5LCBibGFzdCBkdXJhdGlvbiwgYW5kIGJsYXN0IGludGVuc2l0eSAqIGR1cmF0aW9uKSB0aGFuIGJvdGggdGhlIHdhaXRsaXN0IGdyb3VwIGFuZCB0aGUgbWVkaXRhdGlvbiBncm91cC4NCg0KIyMjIyBHcm91cCBEaWZmZXJlbmNlcyBhdCBUaW1lIDMNCg0KT3VyIGNvbnRyYXN0cyBhbmFseXNlcyBhbHNvIHJldmVhbGVkIGdyb3VwIGRpZmZlcmVuY2VzIGF0IFRpbWUgMy4gSG93ZXZlciwgb25seSB0aGUgcmVmbGVjdGlvbiBncm91cCBzaG93ZWQgbGFzdGluZyBwb3NpdGl2ZSBlZmZlY3RzIG9uIGF0dGl0dWRlcyAoc3RpbGwgc21hbGwgZWZmZWN0KSwgYmVoYXZpb3VyYWwgYWdncmVzc2lvbiAoc3RpbGwgc21hbGwgZWZmZWN0KSwgYW5kIGNvbXBhc3Npb24gKHN0aWxsIG1vZGVyYXRlIGVmZmVjdCksIHN1Z2dlc3RpbmcgdGhlc2UgZWZmZWN0cyBhcmUgZHVyYWJsZSBpbiB0aW1lLiBGdXJ0aGVybW9yZSwgdGhlIHJlZmxlY3Rpb24gZ3JvdXAgc2hvd2VkIGEgZGVsYXllZCBvbnNldCBlZmZlY3Qgb24gd2lsbGluZ25lc3MgdG8gaGVscCwgd2hlcmVhcyB0aGV5IHdlcmUgYSBsaXR0bGUgbW9yZSB3aWxsaW5nIHRvIGhlbHAgaW4gdmFyaW91cyBoeXBvdGhldGljYWwgc2NlbmFyaW9zIHRoYW4gdGhlIGNvbnRyb2wgZ3JvdXAuDQoNCiMjIyMgTW9kZXJhdGlvbnMgYXQgVGltZSAyDQoNCkZpcnN0LCBhdHRpdHVkZXMgdG93YXJkIGFnZ3Jlc3Npb24gKE5PQkFHUykgb25seSBtb2RlcmF0ZWQgb25lIHZhcmlhYmxlLCBibGFzdCBkdXJhdGlvbi4gSW4gc2hvcnQsIHdoaWxlIE5PQkFHUyBkb2VzIG5vdCBhZmZlY3QgYmxhc3QgZHVyYXRpb24gaW4gdGhlIGNvbnRyb2wgY29uZGl0aW9uLCBpdCByZWxhdGVzIHRvIGhpZ2hlciBibGFzdCBkdXJhdGlvbiBpbiB0aGUgZGVwbGV0aW9uIGNvbmRpdGlvbi4gQWx0aG91Z2ggdGhpcyByZXN1bHQgaXMgdGhlb3JldGljYWxseSBjb25zaXN0ZW50IHdpdGggdGhlIGxpdGVyYXR1cmUsIGl0IGlzIGxpa2VseSBhIGZhbHNlIHBvc2l0aXZlIGdpdmVuIG91ciBoaWdoIG51bWJlciBvZiB0ZXN0cywgdGhlIGZhY3QgdGhhdCB0aGlzIGlzIHRoZSBvbmx5IHZhcmlhYmxlIHRoYXQgTk9CQUdTIG1vZGVyYXRlcywgYW5kIHRoYXQgdGhlICpwKiB2YWx1ZSBpcyByZWxhdGl2ZWx5IGNsb3NlIHRvIDAuNS4gRnVydGhlcm1vcmUgdGhhdCB2YXJpYWJsZSAoYmxhc3RkdXJhdGlvbiBhbG9uZSkgd2FzIG5vdCBpbiB0aGUgcHJlcmVnaXN0cmF0aW9uLCBzbyB3ZSBtaWdodCBub3QgcmVwb3J0IHRoaXMgZmluZGluZy4NCg0KU2Vjb25kLCBpbXBsaWNpdCBhZ2dyZXNzaW9uIChJQVQpIHNlZW1zIHRvIGhhdmUgbW9kZXJhdGVkIHNldmVyYWwgdmFyaWFibGVzLiBMaWtlIGZvciBOT0JBR1MsIGl0IGFsc28gbW9kZXJhdGVkIGJsYXN0IGR1cmF0aW9uLCBidXQgaW4gYSB0aHJlZS13YXkgaW50ZXJhY3Rpb24gdGhpcyB0aW1lLiBTdXJwcmlzaW5nbHksIGZvciB0aGUgbWVkaXRhdGlvbiBncm91cCwgaW1wbGljaXQgYWdncmVzc2lvbiByZWxhdGVkIHRvIGxvd2VyIGFnZ3Jlc3Npb24sIGJ1dCBvbmx5IHdoZW4gZGVwbGV0ZWQsIHdoZXJlYXMgdGhlcmUgd2FzIG5vIHN1Y2ggaW50ZXJhY3Rpb24gaW4gdGhlIHdhaXRsaXN0IGdyb3VwLiBIb3dldmVyLCBhcyBtZW50aW9uZWQgYmVmb3JlLCB0aGF0IHZhcmlhYmxlIChibGFzdGR1cmF0aW9uIGFsb25lKSB3YXMgbm90IGluIHRoZSBwcmVyZWdpc3RyYXRpb24sIHNvIHdlIG1pZ2h0IG5vdCByZXBvcnQgdGhpcyBmaW5kaW5nLg0KDQpUaGlyZCwgaW1wbGljaXQgYWdncmVzc2lvbiBhbHNvIG1vZGVyYXRlZCBjb21wYXNzaW9uYXRlIGxvdmUsIGFnYWluIGluIGEgdGhyZWUtd2F5IGludGVyYWN0aW9uLiBGb3IgdGhlIHdhaXRsaXN0IGdyb3VwLCB0aGUgZWZmZWN0IG9mIGltcGxpY2l0IGFnZ3Jlc3Npb24gY2xlYXJseSBkZXBlbmRzIG9uIGRlcGxldGlvbjogaW1wbGljaXQgYWdncmVzc2lvbiByZWxhdGVzIHRvIGxvd2VyIGNvbXBhc3Npb24gaW4gdGhlIGNvbnRyb2wgZ3JvdXAgKGV4cGVjdGVkKSwgYnV0IHRvIGhpZ2hlciBjb21wYXNzaW9uIGluIHRoZSBkZXBsZXRpb24gZ3JvdXAgKHVuZXhwZWN0ZWQpLiBIb3dldmVyLCBmb3IgdGhlIG1lZGl0YXRpb24gZ3JvdXAsIHRoZSBlZmZlY3Qgd2FzIGFic2VudCBvciBwYXJ0bHkgcmV2ZXJzZWQuDQoNCkZvdXJ0aCwgaW1wbGljaXQgYWdncmVzc2lvbiBhbHNvIG1vZGVyYXRlZCByZWFjdGlvbiB0aW1lIHRvIHJlbWVtYmVyIGFuIGFsdHJ1aXN0aWMgZXZlbnQsIGFnYWluIGluIGEgdGhyZWUtd2F5IGludGVyYWN0aW9uLiBGb3IgdGhlIG1lZGl0YXRpb24gZ3JvdXAsIGhpZ2hlciBpbXBsaWNpdCBhZ2dyZXNzaW9uIHJlbGF0ZXMgdG8gc2hvcnRlciByZWFjdGlvbiB0aW1lICh1bmV4cGVjdGVkKSwgdW5sZXNzIHRoZXkgYXJlIGRlcGxldGVkLiBIb3dldmVyLCBmb3IgdGhlIHdhaXRsaXN0IGdyb3VwLCB0aGUgZWZmZWN0IHdhcyBhYnNlbnQgb3IgcGFydGx5IHJldmVyc2VkLg0KDQojIyMjIENvbmNsdXNpb24NCg0KSW4gY29uY2x1c2lvbiwgdGhlcmUgc2VlbXMgdG8gYmUgZ3JvdXAgZGlmZmVyZW5jZXMgYXQgVGltZSAyIGFuZCBUaW1lIDMsIGJldHdlZW4gdGhlIGV4cGVyaW1lbnRhbCBjb25kaXRpb25zIGFuZCB0aGUgY29udHJvbCBncm91cC4gSG93ZXZlciwgdGhlIGVmZmVjdHMgaW4gdGhlIHJlZmxlY3Rpb24gZ3JvdXAgYXBwZWFyIG5vdCBvbmx5IHN0cm9uZ2VyLCBidXQgYWxzbyBtb3JlIHJvYnVzdCAoaS5lLC4gdGhleSBhcmUgdGhlIG9ubHkgb25lcyBsYXN0aW5nIGF0IFRpbWUgMykuIEZ1cnRoZXJtb3JlLCB0aGVyZSBhcmUgYWxzbyBzZXZlcmFsIHRocmVlLXdheSBpbnRlcmFjdGlvbnMgYmV0d2VlbiBpbXBsaWNpdCBhdHRpdHVkZXMsIGVnbyBkZXBsZXRpb24sIGFuZCBncm91cCwgYXMgZXhwZWN0ZWQuIFRoZSBuYXR1cmUgb2YgdGhlIGludGVyYWN0aW9ucyBkbyBub3Qgc2VlbSBob3dldmVyIHRvIHBlcmZlY3RseSBhbGlnbiB3aXRoIG91ciBvcmlnaW5hbCBwcmVkaWN0aW9ucy4gQSBkZWVwZXIgZXhwbG9yYXRpb24gb2YgdGhlIG1lYW5pbmcgb2YgdGhlc2UgaW50ZXJhY3Rpb25zIHdpbGwgYmUgcmVxdWlyZWQuDQoNCiMgRnVsbCBDb2RlICYgUmVmZXJlbmNlcyB7LnRhYnNldH0NCg0KVGhlIHBhY2thZ2UgcmVmZXJlbmNlcyBhbmQgdGhlIGZ1bGwgc2NyaXB0IG9mIGV4ZWN1dGl2ZSBjb2RlIGNvbnRhaW5lZCBpbiB0aGlzIGRvY3VtZW50IGlzIHJlcHJvZHVjZWQgaW4gdGhlIHRhYnMgYmVsb3cuDQoNCiMjIC4uLg0KDQojIyBQYWNrYWdlIFJlZmVyZW5jZXMNCg0KYGBge3IgcmVmZXJlbmNlcywgd2FybmluZz1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgcmVzdWx0cz0nYXNpcyd9DQpzZXNzaW9uSW5mbygpICU+JSByZXBvcnQgJT4lIHN1bW1hcnkNCg0KcmVwb3J0OjpjaXRlX3BhY2thZ2VzKHNlc3Npb25JbmZvKCkpDQpgYGANCg0KIyMgRnVsbCBDb2RlDQoNCmBgYHtyIGZ1bGxfY29kZSwgcmVmLmxhYmVsPWtuaXRyOjphbGxfbGFiZWxzKClbIWtuaXRyOjphbGxfbGFiZWxzKCkgJWluJSBrbml0cjo6YWxsX2xhYmVscyhlY2hvID09IEZBTFNFKV0sIGV2YWw9RkFMU0V9DQpgYGA=